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13.  ABSTRACT 


A  study  was  performed  to  analyse  and  minimize  the  effects  of  the  perturbative 
forces  on  a  -24-hour  satellite  orbit  having  an  eccentricity  of  0.6  and  an  inclination 
of  60°  by  selection  of  an  optimum. launch  date.  The  computer  model  used  includes  the 


perturbations  due  to  the  oblate  earth  as  described  by  the  zonal,  tesseral  and 
sectorial  harmonic  series  through  6.  order,  the  solar  and  lunar  third -body  forces,  and 
the  solar-  radiation  pres 3uro  on  a  vehicle  whose  aroa-to-mass  ratio  is  10  ft^/slug,  The 


orbit  was  computed  for  one  year  following  launch  with  orbital  parameter  variations 
presented  as  functions  of  the  launch  date  during  19?0,  The  results  of  tho  study 
Indicate  that  the  total  required  impulse  is  significantly  dependent  on  the  launch  date 
and  oscillatory  behavior  of  the  orbital  parameters  due  to  the  lunar-oyole ,  The  largest 
component  of  the  total  impulse,  due  to  the  plane  change  to  correct  the  nodal  line  and 
inclination,  amounted  to  approximately  66  percent  of  tho  total  yearly  magnitude.  The 
minimum,  total  yearly  Impulse  of  405  m/oec  was  required  when  daily  o  ital  corrections 
were  made  fo.ciuwing  a  30  May  1970  launch,  Uhor.  tho  orbit  was  corrected  at  14-day 


intervals,  in  phase  with  tho  lunar-cycle,  two  local  minimums  wore  found  to  bo  35,5  m/sec ! 
for  a  15  dune  1970  launch  or  360  m/sec  for  a  26  March  19?0  .Launch,  | 
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Preface 

This  report  is  the  result  of  a  study  of  the  effects  of  perturbative 
forces  on  a  specific  high  eccentricity,  high  inclination  24-hour 
satellite.  In  particular,  I  have  attempted  to  determine  the  optimum 
launch  date  in  an  arbitrary' time  frame  of  one  year  on  the  basis  of 
minimum  corrective  impulse  required  to  maintain  the  nominal  orbital 
parameters.  The  work  was  inspired  by  and  in  some  respects  parallels 
the  work  done  by  Arnold  Ricci  (Ref  11)  in  which  a  method  of  making  a 
choice  of  orbital  parameters  to  minimise  the  perturbative  effect  on  the 
ground  track  was  investigated.  It  is  hoped  that  this  study  will  provide 
a  measure  of  value  to  the  future  student  of  perturbative  analysis, 

I  would  like  to  acknowledge  the  special  help  and  guxaance  provided 
by  Prof,  Peter  Bielkowica  who  introduced  me  to  the  fascinating  study  of 
asirodynamios  in  general  and  suggested  this  area  of  special  study. 
Finally,  I  wish  to  express  my  deepest  appreciation  to  my  wife  and  family 
for  their  patience  and  encouragement  during  the  course  of  this  study. 

Allen  L,  Thede 
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'  ...  -  •-  Abstract 

£  study  was  performed  to  analyse  and  minimis©  the  effects  of  the 
perturbative  forces  on  a  24-hour  satellite  orbit  having  an  eccentricity 
of  0,o  and  an  inclination  of:  6-Qr  oy  selection  of  an  optimum  launch  date. 
The  computer  model  used  includes  the  perturbations  due  to  the  oblate 
earth  as  described  by  the  zonal,  tesseral  and  sectorial  harmonic  series 
through  6™  order,  the  solar  and  lunar  third “body  forces,  and  the  soxar 

n 

radiation  pressure  on  a  vehicle  whose  area-to-mass  ratio  is  10  ^/slug. 
The  orbit  was  computed  for  one  year  following  launch  with  orbital 
parameter  variations  presented  as  functions  of  the  launch  date  during 
1970,  The  results  of  the  study  indicate  that  the  total  required  impulse 
is  significantly  dependent  on  the  launch  date  and  oscillatory  behavior 
of  the  orbitax  parameters  due  to  the  lunar-cyle.  The  largest  component 
of  the  total  impulse,  due  to  the  plane  change  to  correct  the  nodal  line 
and  inclination,  amounted  to  approximately  66  percent  of  the  total  yearly 
magnitude.  The  minimum  total  yearly  impulse  of  405  m/sec  was  required 
when  daily  orbital  corrections  were  made  following  a  30  May  1970  launch. 
When  the  orbit  was  corrected  at  14-day  intervals,  in  phase  with  the 
lunar-cycle,  two  local  miniraums  found  to  be  355  m/sec  for  a 
15  June  1970  launch  or  360  m/sec  for  a  26  Larch  1970  launch. 
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MINIMIZATION  OF  PERTURBATIVE  EFFECTS 
ON  A  HIGH  ECCENTRICITY,  HIGHLY  INCLINED 
24-HOUB  SATELLITE  ORBIT 

I*  Introduction 

Synchronous  (24-Hour)  satellites  are  potentially  very  useful  in 
communication  relay  networks,  weather  and  oceanographic  data  collection 
and  ecological  surveys.  Satellites  with  appropriate  orbital  parameters 
can  provide  a  myriad  of  ground  tracks  or  earth  traces  which  are 
characterized  by  figure -8  patterns  (Ref  5)»  These  patterns  can  provide 
optimized  data  collection  or  retrieval  for  opposite  hemispheres  based 
on  the  relative  times  spent  in  various  parts  of  the  ground  track.  The 
potential  application  of  these  qualities  to  various  reconnaissance 
missions  is  also  apparent. 

Satellites  in  orbits  with  large  inclinations  with  respect  to  the 
equatorial  plane  and  having  large  eccentricities  are  subject  to 
significant  perturbations  due  to  the  oblate  earth  and  lunar  and  solar 
attractive  forces.  These  perturbations  will  act,  in  the  long  term,  to 
cause  the  ground  track  to  drift  significantly  from  the  nominal  positions. 
In  order  to  optimize  the  orbital  station  keeping  of  the  satellite, 
consideration  must  be  given  to  the  choice  of  the  optimum  launch  date. 

Equatorial  circular- synchronous  satellite  orbits  have  been 
investigated  extensively  in  the  available  literature.  These  represent 
the  classically  considered  "synchronous"  satellites  which  remain  generally 
over  a  single  point  on  the  earth  surface,  Highly  eccentric  orbits 
(e  >  0,8)  with  orbital  inclinations  of  approximately  30* have  been 
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analysed  by  the  National  'Aeronautics  and  Space  Administration  (NASA)  to 
fulfill  experimenter  requirements  of  solar  aspect  angle  and  position 
relative  to  the  ecliptic  plane  (Ref  12).  Ricci  has  studied  the  effects 
of  general  perturbations  on  24-Hour  satellites  in  order  to  determine  how 
the  injection  parameters  could  be  modified  to  minimize  the  peri ^ bat ions 
in  the  ground  track  (Ref  11). 

This  study  is  concerned  with  the  detailed  analysis  of  a  specific 
orbit  having  an  eccentricity  of  0,6  and  an  inclination  of  60°.  The 
period  of  revolution  is  equal  to  the  sidereal  time  of  the  earth  rotation 
about  its  axis.  The  other  specific  orbital  parameters  were  chosen  to 
provide  a  small  loop  at  the  maximum  southern  latitude,  over  the  Indian 
Ocean.  This  orbit  does  not  correspond  to  any  known  orbit  at  the  present 
time. 

The  purpose  of  the  analysis  is  to  determine  the  optimum  launch 
date  during  the  arbitrary  calendar  year  1970,  The  optimum  date  will  be 
chosen  on  the  basis  of  minimum  corrective  impulse  to  provide  nominal 
station  keeping  for  a  period  of  one  year,  Tho  computer  model  used  in 
the  analysis  applies  the  Enoke  technique  of  integration  and  was  written 
to  be  used  on  the  CDC  6600  computer  at  the  Wright -Patterson  AFB 
Computer  Selena.  Center, 

The  oblate  earth  has  been  modeled  by  using  the  J„  through  J,  terms 

Z  0 

of  the  zonal  harmonics  and  includes  the  tesseral  and  seotorial  harmonios 
through  the  6  order  terms.  The  effects  of  tho  lunar  and  solar  third- 
body  attractive  forces  as  well  as  the  radiation  pressure  due  to  solar 
photons  have  been  inoluded  in  the  model.  The  vehicle  is  assumed  to  be  a 
spherically  symmetric  perfectly  absorptive  body  with  an  area-to-mass 
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ratio  of  10  ft  /slug.  The  perturbative  forces  due  to  atmospheric 
drag  on  the  vehicle  and  the  third-body  attractions  of  the  planets  have 
been  neglected  since  their  magnitudes  are  insignificant  compared  to  the 
other  sources  considered. 

Several  other  basic  assumptions  have  been  made  during  the  course 
of  this  study.  It  is  assumed  that  the  gravitational  field  of  the  earth 
is  accurately  modeled  by  the  harmonic  series  used.  It  is  assumed  that 
the  satellite  has  been  injected  at  perigee  of  the  initial  orbit  from 
some  available  launch  facility.  The  study  does  net  include  analysis 
of  the  injection  phase.  It  is  assumed  that  the  satellite  attitude 
control  system  operation  does  not  significantly  perturb  the  orbit. 

The  analysis  and  study  results  are  presented  in  the  chapters 
which  follow.  The  general  analysis  approach  is  discussed  in  Chapter  II 
with  detailed  formulations  developed  and  presented  in  the  appended 

i 

sections.  The  results  of  the  study  phases  are  presented  largely  in 
graphical  form  in  Chapter  III,  The  final  chapter  provides  a  summary 
of  conclusions  and  makes  recommendations  for  further  study. 
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XI.  Analysis  Approach 

This  chapter  of  the  study  presents  the  methodology  used  during  the 
analysis  beginning  with  a  definition  of  the  general  approach.  The  second 
section  is  devoted  to  the  discussion  of  integral  parts  of  the  computer 
mode],  used,  The  aspects  of  choosing  initial  conditions  for  the  various 
launch  dates  are  discussed  in  the  third  section.  Finally,  the  last 
section  presents  the  technique  used  to  determine  tnc  impulse  necessary  to 
provide  station  keeping. 

General  Discussion 

The.  approach  to  the  analysis  of  the  orbit  consisted  of  three  specific 
phases.  The  first  phase  involved  the  study  of  daily  variations  in  the 
orbital  parameters.  The  variations  were  calculated  as  functions  of  each 
of  the  perturbation  sources  for  a  typical  orbital  revolution  and  for  the 
total  perturbative  force.  The  second  phase  included  the  calculation  of 
the  orbit  for  a  full  year  following  injection  at  perigee  on  tho  launch 
date.  Successive  launch  dates  were  chosen  throughout  the  calendar  year 
of  19?0,  The  final  phase  provided  for  the  select ion  of  an  optimum  launch 
date  by  calculation  of  the  impulse  required  fur  station  keeping  given 
various  launch  dates, 

The  specific  orbit  which  was  chosen  for  analysis  las  the  following 
orbital  parameters) 

a  -  6,610725  R_ 

O 

e  •  0,6 

i  -  60° 

u  *  135° 

O  *  0°  (for  0h  l  Jan  19?0  launch  time) 
TP  -  23h  56®  4.09054® 
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The  choice  of  the  longitude  of  ascending  node  ft  for  launch  dates  other 
than  Tq  *»  1  Jan  1970  are  discussed  in  the  third  section  of  this  chapter 
■on  page  15,  with  the  initial  launch  conditions  illustrated  in  Fig,  7. 

Moat  of  me  orbital  calculations  in  the  anal yn Is  wex'e  performed  in  the 
geocentric  equatorial  frame  of  reference  system  XYZ,  defined  by  the  unit 

A  \  A 

vectors  I,  J  and  K,  The  perifocal  system  of  coordinates  PQW  was  used 
to 'facilitate  calculations  in  the  orbital  plane,  The  geometry  of  thes 
frames  of  reference  is  shown  in  Fig,  1, 

Computer  Model 

fern  nee  Orbit.  The  reference  orbit  was  calculated  according  to 
the  development  in  Appendix  B,  The  orbital  parameters  used  for  the  first 
reference  orbit  are  the  nominal  values  given  previously.  The  first  step 
in  the  calculation  is  to  determine  the  values  of  the  eccentric  anomaly 
E  by  solution  of  the  Kepler  equation 

mw  -  (E  -  0  sin  E)  (1) 


where  m  is  integer  step  and  w  is  the  desired  time  interval  between 
equally  apaoed  data  points  around  the  orbit.  An  interval  of  w  ■  0,1  hr 
was  used  la /.this  analysis. 


^  position  and  velocity  components  of  the  reference  orbit  were 
then  calculated  by  using  the  general  vector  equations 


r  • 

v  • 


A  A 

r  cca  v  P  +  r  siny  Q 


ein  v  P  +  (e  +  cos  u )  5]  (3) 


Theoe  equations  are  transferred  from  the  perifocal  reference  system  HJV 
to  the  equatorial  coordinate  system  XX2  in  Appendix  B, 
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Bert ur bat Ions. 

l)  Oblate  Earth.  The  components  of  the  perturbative  accelerations 
due  to  the  oblate  earth  were  calculated  using  the  form  of  the  gravitational 
potential  expressed  as 


0 


K 


_e 

r 


6 

-El 

n-2 


yr/^(sin  «)' 

.k-0 


[enkoos(JtA.E)  +  E^slnfkXj,)] 


W 


A  6^  coder  series  is  used  to  assure  sufficient  aoouracy  while  not  being 
unwieldy  in  computer  calculations.  The  gradient  of  0  can  be  decomposed, 
as  in  Appendix  A,  to  yield  the  zonal,  tessoral  and  sectorial  harmonic 
functions  which  desoribe  the  perturbable  contributions  to  the  gravitational 
acceleration  of  a  body  in  earth  space.  These  functions  ace  given  by  Eqs 
(22)  through  (24)  in  Appendix  A,  end  Illustrated  in  Figs.  2  through  4 

■W 

for  the  first  orbit  following  Injection  at  perigee  at  0  1  Jan  1970, 

It  can  be  seen  that  the  tonal  contributions  are  tho  most  significant 

with  peak  perturbations  near  perigee  being  on  the  <«r  of  lO*^  B#/W2, 

The  second  tost  significant  components  are  those  of  the  sectorial  harmonic* 

which  describe  the  equatorial  elllptieity,  Beak  sectorial  contribution* 

are  on  the  order  of  10”^  a  /hr2,  The  tesseral  harmonics,  which  further 

e 

describe  the  nature  of  the  go  old,  have  peaks  on  the  order  of  10*^  ft#/hr2. 

The  tesseral  )oxconics  contribute  very  little  to  the  variation  of  orbital 
parameters,  ts  will  be  shown. 
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Z)  Third-Body  Attraction.  The  contributions  to  the  perturb¬ 
ative  accelerations  due  to  third-body  attractive  forces  of  the  &000  and 
sun,  in  the  geocentric  equatorial  system,  were  calculated  by  the  general 
form  of  Eqs  (69)  and.  (?0)  given  in  Appendix  Aj 


•c 


The  contributions  for  the  orbit  beginning  at  0^  1  Jan  1970  due  to  the 
lunar  and  oolax  attractions  are  shown  in  Figs,  5  and  6  respectively. 


Although  the  maximum  accelerations  axe  on  the  carder  of  10  and 
10  0  K  /hr  ,  they  are  present  for  longer  periods  of  tins  during  the 
orbit.  Their  duration  causes  perturbative  displacements  in  the  orbit 
which  exceed,  the  effects  of  the  sectorial  harr.onics  and  approach  those 
due  to  even  the  .:oaal  harmonies. 


3)  Solar  Bad .iat ion  Presaure .  Finally ,  Use  contribution  due 
to  tlte  solar  radiation  is  calculated  by  the  equation 


t  v 
s  40 


A 

R 
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Although  the  contributions  duo  to  this  perturbative  source  axe  not  shown 

n  * 

they  era  on  the  cud  or  of  10'  R  /hr*  fc«.  an  area-to-aasn  ratio  >-  - 

I?  3 

10  ft  /slug,  Since,  in  the  heliocentric  system,  the  satellite  is  at  & 
relatively  constant  distance  frees  the  sun,  the  perturbative  acceleration 

L 

is  also  nearly  constant,  The  x,  y,  and  c -compose  his  for  the  0  1  Jan  X9?0 
orbit  are  approximately  -1.15  *  10~' ,  5,6?  x  10  ^  and  2.^5  x  10  1  R^/hr2 
respectively , 


n 


8  10  12  16 


20  22  2-* 


Hour  in  Orbit 


olar  third -&ody  Fort ur bat Sva  Acceleration  Co»pcbentc 
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The  value  of  10  ft  /slug  is  considered  to  be  a  maximum  value  for 

~  for  most  typical  satellites,  x.:  ^act,  this  value  is  higher  than 

envisioned  for  many  vehicles  designed  with  large  solar  arrays  necessary 

for  some  ion-propulsion  applications.  The  factor  f  can  ba  considered 

3 

A 

to  be  identically  equal  to  unity  for  this  ~  when  the  vehicle  is  in  the 
sunlit  portion  of  the  orbit,  and  e^  «tl  to  v-ero  when  it  enters  the 
shadowed  areas  behind  the  earth. 


4)  Neglected  Perturbations,  The  most  significant  perturbation 
which  has  teen  neglected  in  this  analysis  is  that  due  to  the  atmospheric 
drag  on  the  vehicle,  particularly  at  perigee,  It  can  be  shown,  however, 
that  its  contribution  is  negligible  compared  to  those  already  discussed. 
The  total  Gecelerative  perturbation  along  the  path  of  the  orbit  is  given 
by  the  general  relation 


(7) 


The  greatest  drag  force  occurs  at  perigee  where  the  nominal  vexocity  v 

is  6,149  Km /sec  *  20,175  ft /sec  and  the  altitude  is  10488  km  » 

£ 

34,408  x  10  ft.  The  d'nsity  o  at  this  altitude  is  certainly  no 
greater  than  lO”1^  slug/ft"’  (Ref  7*200),  Considering  the  case  of 

i  o 

-  »  10  ft  /slug  and  assuming  that  the  drag  coefficient  maximum  is 
m 

CL  -  2,  the  deceleration  r  is  no  greater  than  2.52  x  10~10  R/hr2, 

D  6 

This  value  as  at  least  three  orders  of  magnitude  loss  than  the  tesseral 
component  which  if  the  smallest  contribution  considered  in  the  analysis. 

In  general,  a  less  significant  perturbation  results  from  the  third- 
body  attraction  of  near-by  planets,  specif ioaliy  the  larger  ones.  The 
approximate  values  of  these  perturbative  accelerations  can  be  obtained 
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oy  using  a  simplified  form  of  Eq  (5)  gi.en  as 


‘Pi 


-  K 


max 


Pi  r  3 
Pi 


(8) 


tfhsre  is  the  gravitational  parameter  of  the  planet  of 

mass  The  maximum  value  of  the  perturbative  force  will  be  reached 

for  a  maximum  redial  distance  of  the  vehicle  from  the  earth  r  and 

max 

a  minimum  distance  between  the  vehicle  and  the  perturbing  planet,  r^. 
Assuming  the  planets,  including  earth,  are  in  the  equatorial  plane,  the 
minimum  value  of  r  ^  is  at  closest  approach,  or  approximately  the  diff¬ 
erence  between  the  semi-major  axes  of  the  orbits  of  the  earth  and  planet 
about  tha  sun. 


The  maximum  third-body  perturbation  due  to  the  nearest  planets  based 
on  closest  possible  approach  are  tabulated  in  Table  1  below.  The 
planetary  masses,  semi-major  axes  and  the  minimum  approach  distances  axe 
also  provided.  The  perturbations  are  based  on  a  maximum  vehicle  distance 
of  r  -  10.577  R  . 

max  e 


Table  I,  Third -Body  Perturbations  of  Selected  Planets 


Planet 

. . 

Mass 

c,) 

Semi-major 
Axis  (a,u, ) 

rpl 
(a.u, ) 

•  V 

rmax„ 

(Re/hr^) 

Venus 

0.61485 

0,723332 

0,276668 

627.94  xl0~12 

Kars 

0,10770 

1.523691 

0,523691 

12.24  xl0~12 

Jupiter 

317.890 

5.202694 

4,202694 

69.89  xl0~12 

Saturn 

95.1200 

9.538836 

6,538836 

2,49  xl0~12 
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Launch  Date  Variation 

.In  order  to  maintain  the  same  relative  conditions  of  solar  aspect 

for  the  satellite  with  respect  to  the  desired  earth  trace,  it  is  necessary 

to  adjust  the  longitude  of  ascending  node  for  each  launch  date,  Tha 

adjustment  is  based  on  the  launch  (injection  at  perigee)  at  0^  on  the 

date  of  interest.  Since  the  value  at  0^  1  Jan  19?0  is  »  0®,  the 

o 

value  for  any  other  time  of  launch,,  as  shown  in  Pig,  7,  is 


i 


(9) 


where  6  is  the  Greenwich  sidereal  time  at  0n  1  Jan  1970,  which  is 
6h  40m  55.338s  «  100.25°,  and  6gi  is  the  Greenwich  sidereal  time 
at  0h  on  the  launch  date  (fief  ltlO),  The  dates  shown  in  Fig,  7 
represent  launch  dates  for  initial  orbital  positions  shown  in  the 
heliocentric  equatorial  system  X '  Y '  Z '  , 


The  longitude  and  latitude  of  the  vehicle  with  respect  to  the  earth 
is  given  by  the  relations 

AE  “  01  "  egi  “  we  “  V  (10) 

0  -  sin"1  (|)  (11) 

where  a  »  is  the  right  ascension  of  the  vehicle  and  (t-tQ) 

Is  the  time  from  perigee.  Given  the  nominal  orbital  parameters,  Eqs  (10) 
and  (ll)  can  be  i  od  to  generate  the  ground  traok  or  earth  trace  as 
shown  in  Fig,  8,  The  points  along  the  ground  track  in  Fig,  8  represent 
intervals  of  one  hour,  which  demonstrates  the  utility  of  an  orbit  of  this 
typo  for  collection  of  data  in  the  general  area  of  the  small  loop. 
Approximately  one  third  of  the  total  orbital  time  is  spent  in  the  general 
area  of  this  loop. 
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Following  the  orbital  calculations  for  each  of  the  selected  launch 
dates  using  the  variable  longitude  of  ascending  node,  the  perturbed  ground 
track  can  be  calculated  using  Eqs  (10)  and  (ll),  The  value  of  0^, 
however,  must  be  determined  specifically  for  the  time  of  perigee  passage 
for  the  orbit  of  interest. 


Orbital  Correction 

Following  the  necessary  analysis  of  orbital  parameter  variation  as 
a  function  of  launch  date,  a  series  of  orbital  integrations  were  performed 
to  determine  the  optimum  launch  date.  Detailed  analysis  was  made  only 
for  the  selected  general  area  in  order  to  minimize  the  computer  time 
required.  The  method  used  to  correct  the  orbit  is  discussed  in  more 
detail  in  Appendix  B  although  basically  it  involved  making  four  impulsive 
thrusts  during  each  orbit  being  corrected. 

The  first  correction  was  an  orthogonal  impulse  made  near  the  lowest 
point  L  in  the  orbit  below  the  equatorial  plane  to  correct  the  values 
of  inclination  and  longitude  of  ascending  node.  This  correction  is 

A  i  vT  An  vT  sin  iT 

Av  -  - k -  -  - - k - a -  (12) 

cos  (^  -  ifc)  ,  sin  (  ) 

-l  /An 

where  ^  "  tan  I  —  sin 


The  second  impulse  was  made  at  apogee  to  correot  the  apsidal  line 
or  argument  of  perigee.  It  was  made  normal  to  the  velocity  vector, 
toward  the  center  of  attraction  and  has  a  magnitude  given  by 


Ay 


\  6a  va 


°  2  aA  +  rA  OOB  (»A  -  »#) 


(13) 
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The  third  impulse  was  also  made  at  apogee  but  in  a  direction  which 
is  tangential  to  the  orbit.  This  correction  provides  a  partial  correction 
of  the  eccentricity  and  semi-major  axis  and  is  based  on  restoring  the 
perigee  distance  to  the  nominal  value.  This  impulse  is  given  by 


Av 


ta 


AaAKe 
2  aA2  VA 


(14) 


The  final  tangential  impulse  was  given  at  perigee  to  fully  restore 
the  nominal  eccentricity  and  semi-major  axis  and  is  given  by 


Av 


tp 


Ae  V* 

_ E _ E_ 


3  (e^  +  cos  ip 


(15) 


where  the  primed  values  are  those  just  prior  to  the  impulse  being  made. 

The  choice  of  the  optimum  launch  date  finally  depends  on  the 
criteria  of  total  impulse  necessary  for  orbital  correction.  Although 
the  ultimate  design  of  the  thruster  assembly  to  perform  the  correction 
is  l  yond  the  soope  of  this  study,  an  understanding  of  the  basic 
requirements  is  of  interest.  Due  to  the  oscillatory  nature  of  the 
parameter  changes,  the  impulse  required  for  daily  correction  will  be 
somewhat  higher  than  if  the  correction  were  made  periodically  as  a  function 
of  the  lunar  cyole.  The  final  choice  of  the  corrective  technique  depends 
on  the  acouracy  which  is  required  for  station  keeping, 
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III,  Results 

The  results  of  the  various  phases  of  this  study  are  discussed  in 
this  chapter  with  the  aid  of  graphical  presentations  of  the  data  acquired. 
The  effects  of  the  perturbations  on  the  orbital  parameters  during  a  typical 
orbit  are  discussed  in  the  first  section.  The  longer  term  effects  of 
these  perturbative  sources  are  presented  in  the  second  section  as 
functions  of  the  date  in  orbit  and  launch  date.  Finally,  the  optimum 
launch  date  is  selected  on  the  basis  of  minimum  corrective  impulse  being 
required  to  maintain  nominal  orbital  parameters. 

Hourly  Variation  of  Orbital  Parameters 

The  effects  of  the  various  mechanisms  of  perturbation  on  the  orbital 
parameters  are  shown  for  the  first  orbit  after  injection  at  0^  1  Jan  1970 
in  Figs,  9  through  13 ,  The  individual  curves  for  the  various  components 
of  perturbation  are  explicitly  the  effeot  due  to  those  components  alone, 
while  the  total  effect  is  due  to  the  sum  of  the  accelerations  considered 
together.  It  can  be  seen,  however,  that  an  accurate  approximation  can  bo 
determined  by  summing  the  effects  due  to  each  individual  contribution, 
if  a  component  suoh  as  solar  radiation  pressure  or  tesseral  harmonic  is 
not  shown,  there  is  a  negligible  change  in  the  nominal  value  of  the 
parameter  considered. 

Figure  9  illustrates  the  effect  of  the  various  pert ur hat ions  on  the 
semi-major  axis  during  this  typical  orbit,  Hots  thar  the  moat  significant 
effeot  is  caused  by  the  zonal  harmonics.  The  sharp  s  piker  occurring  near 
perigee  will  also  be  seen  in  the  variation  of  other  orbital  parameters. 

As  discussed  by  Shuts  (Ref  12«15),  the  effects  of  these  spikes  will  bo 
averaged  out  following  perigee  passage.  The  variations  in  the  semi-major 
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axis  due  to  the  zonal  harmonics  Kill  then  he  of  the  same  order  of  ' 
magnitude  as  compared  to  the  other  contributions  at  the  perigee  point. 

The  effects  of  the  earth  oblateness  will  remain  essentially  as  shown  for 
any  orbit  given  any  launch  date,  since  the  orbit  is  constant  relative 
to  the  surface  of  the  earth.  This  will  also  be  true  for  the  other  orbital 
parameters  discussed. 

The  next  most  significant  contributions  which  cause  variation  in 
t he  semi-major  axis  are  the  lunar  and  solar  third-body  attractions. 

Unlike  the  oblateness  effects,  these  variations  Kill  depend  on  the 
positions  of  the  moon  and  the  sun  with  respect  to  the  earth.  The  variations 
shown  in  Fig,  9  are  therefore  only  typical  in  magnitude.  The  variations 
due  to  solar  radiation  pressure , with  ~  •»  10  1  /slug,  and  the  tessewal 
harmonics  (not  shown)  are  seen  to  be  minimal  in  comparison  with  the  other 
contributions. 

The  variation  of  the  eccentricity  is  illustrated  in  Fig,  10  where 
it  is  no4  od  that  the  lunar  attraction  is  the  most  significant  perturbation 
source.  The  solar  attraction  offsets  this  effect  somewhat  far  this 
typical  orbit,  as  do  the  zonal  harmonics.  The  effects  due  to  the 
tesseral  and  sectorial  harmonics  and  the  solar  radiation  pressure  are 
eigaif  leant  compered  to  those  due  to  the  sources  shown. 
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The  inclination  is  also  affected  Boat  significantly  by  the  lunar, 
solar  and  aonal  oontri^tions  as  shown  in  Fig «'  11.  the  sectorial  harmonica 
and  solar  radiation  pressure  cause  less  significant  variation  although 
cannot  be  neglected,  the  solar  attraction  for  the  orbital  case  shown 
causes  an  increase  in  the  inclination  which  is  twice  the  Increase  caused 
by  the  canal  iuxsonics,  while  the  lunar  attraction  cancels  scat  of  its 
effect. 

the  argues nt  of  perigee  u>  is  affected  aost  significantly  by  the 
stoiial  haraonles  as  shown  in  Fig*  12*  The  solar  and  lunar  attractive 
influences  tend  to  cancel  for  this  particular  case,  but  as  noted  earlier, 
will  act  generally. 
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Hour  la  Orbit 

Figure  13,  Variation  of  Longitude  of  Ascending  Node  with  Hour  in  Orbit 


finally,  Fig,  13  above  illustrates  the  hourly  change  in  the 
longitude  of  the  ascending  node  for  the  orbit  following  injection  at  0h 
1  dan  1970  where  £«ft  *  O'* ,  The  $ca*J,  solar  and  lunar  perturbations  are 
seen  to  be  the  scat  effective  in  the  variation  of  this  parausetcr. 


&&  MMlSB 

I'mnrMns*  se  E&to  in  Orbit,  the  orbit ,  uas  calculated  for 
a  period  of  or ,  /oar  following  the  e*iected  launch  dates  which  mm 
chosen  during  the  period  of  1  Jan  1970  through  1  dan  1971.  The  results 
of  selected  orbital  coagulations  which  are  illustrated  in  Figs,  14  through 
19  represent  the  changes  of  the  orbital  jaraaoters  as  a  function  of  the 


tins  in  orbit*  The  separate  curves  shown  in  each  figure  represent 
different  launch  dates  Te  (injection  of  vehicle  into  orbit  at  perigee). 
Ifee  specific  launch  date  oases  shown  were  eelecUf  l?  illustrate  in 
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general  the  extremes  of  the  parameter  variation  which  result  regardless 
of  launch  date.  The  individual  curves  are  based  on  the  value  of  each 
parameter  at  successive  perigee  points  of  woe  orbit. 

The  short-period  oscillatory  behavior  evident  in  these  figures  is 
the  effect  of  the  lunar-cycle,  while  the  longer  period  oscillations  are 
due  to  the  solar-cycle.  These  oscill ~ A * ons  are  superimposed  on  the 
general  variation  trend  which  Is  due  to  the  oblateness  of  the  earth. 

Figure  14  illustrates  the  variation  of  the  semi-major  axis  as  a 
function  of  tins  in  orbit  following  injeotion  on  the  dates  indicated. 

The  daily  orbital  variation  as  shown  for  T  •  2?4  includes  the  effec  jf 
the  lunar  cycle,  while  the  cases  for  T  •*  69  and  Tq  «  182  illustrate  the 
longer  period  variation  as  a  result  of  smoothing  the  curves,  The  three 
cases  shown  represent  an  envelope  into  which  most  other  launch  oases 
fall.  The  smoothed  carves  wort  used  here  to  avoid  a  cluttered  represent¬ 
ation  and  will  act  be  used  for  other  parameters. 

The  variation  of  eccentricity  is  shv,*a  in  Fig.  15  for  several  launch 
dates.  The  least  variation  in  the  nominal  "&lue  is  obtained  ter  a  laustch 
date  of  T0  *  69  and  the  maximum  change  is  for  a  launch  at  Tq  <*  IBS.  The 
general  trend  in  decreasing  the  eccentricity  causes  the  orbit  to  circul¬ 
arise  in  Vise, 

Figure  16  illustrates  the  varieties  of  inclination  for  a  selected 

group  of  launch  dates.  The  greatest  variation  is  obtained  for  a  launch 

date  of  T  *  51.  while  for  T  *305,  the  inclination  oscillates  about  the 
p  *  *  c 

iU»inal  value  of  60°,  In  general,  the  perturbs! iciis  total  to  increase 
the  orbital  inclination. 
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The  general  trend  of  the  perturbations  lead  to  an  increased  value 
of  the  argument  of  the  perigee  as  shown  in  Fig  1?,  The  greatest  decrease 
in  u  is  obtained  for  a  launch  date  of  Tq  -  121  while  the  maximum  Increase 
is  reached  for  Tq  *=  274,  The  variation  of  the  apsid&l  line  is  most 
dependent  on  the  variations  in  the  orbital  inclination  and  eccentricity. 


■M 


Since  the  initial  values  of  the  longitude  of  ascending  node  ft  varies 
with  the  launch  date,  Fig,  18  represents  the  variation  of  Aft  ,  In 
general,  the  trend  of  westward  regression  is  noted  with  a  minimum  drift 
resulting  for  Tq  «  182  and  maximum  drift  at  Tq  «  51,  The  ..ariation  due 
to  the  lunar-cycle  is  absent  from  these  curves  and  the  solar*- cyclo  is 
only  vaguely  apparent. 


A*** 
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The  variation  of  the  perigee  distance  is  shown  in  Fig,  19  and 

illustrates  the  close  dependence  on  the  variations  in  the  eccentricity 

as  shown  in  Fig,  15,  Note  that  the  minimum  change  in  r  is  obtained 

P 

for  T  “  69  while  the  maximum  Variation  is  reached  at  T  182, 

0  o 


Parameter  Dependence  on  Launch  Date.  In  order  to  more  dearly 
determine  the  effect  of  the  launch  date  choice  on  the  variation  of  the 


orbital  parameters,  the  data  has  been  presented  as  a  function  of  launch 


...  • 

*■  x.r 


date  in  Figs,  20  tiirough  24,  The  separate  curves  in  each  figure 
represent  the  variations  after  90,  180,  270  and  36O  days  in  orbit.  These 


plots  dearly  indicate  the  launch  dates  which  will  provide  the  minimum 
changes  in  given  orbital  parameters.  These  characteristics  In  themselves, 
however,  cannot  be  directly  used  in  selecting  the  optimum  launoh  date 
as  will  bo  shown  in  the  following  section. 
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Since  the  serai-major  axis  varies  in  a  very  small  envelope  regardless 
of  the  launch  date,  its  variation  as  a  function  of  launch  date  is  not 
shown  in  this  section.  Figure  20  represents  the  variation  of  eccentricity 
as  a  function  of  launch  date  which  clearly  Indicates  a  maximum  rate  of 
change  at  a  launch  date  of  Tq  »  182.  There  are  two  launch  dates  which 
provide  minimum  changes  and  .are  at  Tq  «  60  and  To  -  30Q,  The  minimum 
changes  represent  a  decrease  in  the  eccentrioity  by  2.42  percent  while 
the  .maximum  decrease  represents  a  7,85  percent  change. 

The  variation  of  the  inclination  is  shown  in  Fig.  21  which  illustrates 
that  there  are  two  launch  dates  which  provide  minimum  changes  during  the 
year  in  orbit  while  there  are  also  two  dates  which  provide  maximum 
changes,  The  maximum  at  To  »  60  represents  an  increase  in  the  inclination 
hy  3.5  percent  ,  while  the  maximum  at  Tq  -  222  is  an  increase  of  only 
0,83  percent.  The  minimum  at  T0  *  165  is  only  a  0.45  percent  increase 
while  that  at  Tq  represents  a  negligible  change  by  oscillating  about  the 
nominal  value  of  60°, 

The  variation  in  tta  argument  of  perigee  shown  in  Fig.  22  indicates 
two  possible  launch  dates  which  provide  essentially  no  ohanget  TQ  *  70 
and  T  -  175 .  A  maximum  increase  by  3,3  percent  is  found  f or  T  *  280 

w  O 

while  the  maximum  decrease  of  1,48  percent  is  reached  at  TQ  •  125, 

Figure  23  illustrates  the  variation  in  the  longitude  of  the  ascending 

node.  The  greatest  magnitude  is  a  2,3  percent  decrease  at  T  *  50 

o 

(percentage  based  on  a  nominal  value  of  360°),  The  minimum  change  is  a 
1,61  percent  decrease  at  -  176. 
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Tidily,  the  variation  in  the  perigee  distance  is  shown  in  Fig.  24, 
The  maximum  increase  by  11,55  percent  at  Tq  "  182  represents  a  very 
significant  change.  The  minimum  changes  at  Tq  -  70  and  Tq  *  300 

represent  3.62  percent  increases  in  the  nominal  value  of  r  , 

P 

Earth  Trace  Variation.  The  effect  of  the  orbital  parameter  changes 
on  the  earth  trace  after  one  year  in  orbit  is  shown  in  Fig,  25,  The  two 
cases  shown  were  selected  to  illustrate  two  extreme  examples  of  drift. 
The  trace  for  TQ  *  121  represents  the  maximum  westward  drift  of  the 
loop  In  the  southern  hemisphere,  while  the  case  for  Tq  *  274  represents 
the  maximum  eastward  drift. 

The  relative  size  and  position  of  the  small  loop  can  be  roughly 
determined  by  referring  to  Fig,  26,  which  illustrates  the  longitudes 
of  the  most  easterly  point  and  crossover  point  as  functions  of  the 
launch  dates.  Figure  2?  illustrates  the  latitudes  of  the  nest  easterly 
end  crossover  points  in  a  similar  manner. 


Figure  35.  C'Gfijwfcr l«a  of  Barth  Trace*  of  Perturbed  and  I’caiaal  Orbit*  After  One  Tear 


ta.%1  tv'Xm  {&*%)  t*  Longitude  (deg) 


GA/KC/73-7 


vxf.  26 .  Longitudes  of  S&&11  Loop  as  Fu&otioa  of  LauacL  Date 
After  One  Year 


Figure  27*  Latitude*  of  Sisal i  Loop  m  Function  of  iauoen  Lata 
After  Coe  Year 
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Orbital  Correction 

The  total  impulse-  required  for  dally  correction  has  been  compared 
with  total  impulse  required  for  correction  every  14  orbits,  Figure  28 
illustrates  the  magnitudes  of  impulses  required  for  the  two  options  studied. 
As  Indicated  in  the  figure,  the  area  of  interest  is  the  range  of  launch  dates 
providing  for  minimum  corrective  impulse.  Except  where  short-term 
oscillations  result  in  continually  increasing  (or  decreasing)  orbital 
parameters,  the  daily  corrective  impulses  are  greater  in  magnitude  than 
those  made  periodically. 

The  largest  component  of  the  total  corrective  impulse  required  is 
due  to  the  orthogonal  impulse  Av  required  to  correct  the  line  of  nodes 

rf 

and  inclination.  The  magnitude  of  Avf  is  most  heavily  uependent  on  the 
change  required  in  the  longitude  of  ascending  node  rather  than  the 
inclination.  The  minimum  value  for  this  impulsive  component  is  obtained 
at  Tq  ■  1?0  with  268  m/sec  and  244  m/sec  being  required  for  the  daily  and 
the  14-day  interval  corrections  respectively.  These  values  are  nearly 
66  percent  of  the  total  impulse  required  for  correction  after  a  year  in 
orbit. 

The  second  most  significant  contribution  to  the  total  yearly  impulse 
is  "onerally  due  to  the  normal  impulse  necessary  to  correot  the  argument 
of  perigee  or  &p3idal  line.  The  minimum  value  of  Av^  for  daily  correc¬ 
tion  is  at  50  m/sec  at  »  165.  If  corrections  are  made  at  14-day 
intervals,  however,  two  local  minimums  are  found.  The  minimum  at  T  »  I65 

v 

is  only  20  m/V°*c  while  the  one  at  *  *  00  is  40  m/sec.  The  normal  impulse 
approaches  the  magnitude  of  the  orthogonal  impulse  for  a  launch  date  of 
T  »  2?4  when  it  reaches  a  maximum  value  of  approximately  200  m/sec, 


Components  of  Corrective  Impulse  Required  to  Maintain  Nominal  Orbital  Parameters 
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As  is  expected,  the  shape  of  the  curve  depicting  Av^  is  dependent 
primarily  on  the  absolute  magnitude  of  the  variation  of  the  argument  of 
perigee. 

The  tangential  impulses  Av^  and  Av^  are  generally  the  smallest 
components  of  the  total  yearly  correction  and  are  primarily  dependent  on 
the  magnitude  of  the  variation  in  the  eccentricity.  The  impulse  made  at 
apogee  is  necessarily  largest  of  the  two  since  it  depends  inversely  on 
the  orbital  velocity.  The  apogee  impulses  are  four  to  five  times  larger 
than  those  made  at  perigee  and  are  directed  in  the  opposite  direction 
relative  to  vehicle  motion.  The  minimum  values  for  the  tangential  impulses 
are  obtained  for  »  60,  The  magnitudes  of  Av^a  for  daily  and  14-day 
interval  corrections  are  40  m/sec  and  25  m/sec  respectively.  The  minimum 
value  of  Av+p  is  6  m/sec  for  either  technique  of  correction. 

The  minimum  total  impulse  required  to  maintain  the  nominal  orbit 
using  daily  corrections  is  obtained  for  a  launch  date  of  Tq  «  150  and  has 
a  magnitude  of  40 5  m/sec.  The  required  impulse  is  less  than  410  m/sec  for 
a  range  of  launch  dates  from  Tq  *  130  to  Tq  «*  165,  Two  local  minimums 
are  obtained  if  corrections  are  applied  at  14-day  intervals.  These 
points  are  reached  for  To  «  85  and  T  a  166  and  indicate  minimum  required 
impulses  of  36O  m/sec  and  355  m/sec  respectively.  The  technique  of  using 
14-day  intervals  for  corrections  provides  reductions  in  the  total  impulse 
of  14,3  percent  and  13,4  percent  for  Tq  »  85  and  T0  «  166  respectively 
as  compared  to  those  for  daily  corrections. 

The  foregoing  analysis  of  the  minimum  impulse  required  to  maintain 
the  nominal  orbit  should  be  contrasted  with  that  required  if  the  worst 
date  were  chosen  for  launch  during  the  calendar  year  of  1970,  The  total 
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impulse  required  for  a.  launch  on  Tq  «  300  would  be  approximately  558  m/sec 
if  the  orbit  were  corrected  on  a  daily  basis.  The  impulse  would  be  nearly 
this  value  if  corrected  at  14-day  intervals  since  the  orbital  parameters 
change  rapidly  and  the  lunar  oscillations  are  not  as  effective.  This 
value  is  49,4  percent  higher  than  for  the  daily  corrected  optimum  case 
and  55*9  percent  greater  than  the  14-day  interval  case* 


body  perturbations  are  second  only  to 


harmonics.  The  radiation 


axea-to-raass  ratio  of  10  ft  /slug  is  generally  more 


pressure  ior  tne 
significant  than  the  sector! 
probably  be  neglected  for  this  type  of  study. 


rmonics . 


he  tesseral  harmonics  could 


The  study  has  resulted  in  the  determination  of  optimum  launch  dates 
to  minimize  the  total  station  keeping  impulse  required  for  one  year  in 
orbit.  The  optimum  dates,  however,  were  found  to  depend  on  the  specific 
technique  of  impulse  application,  or  more  correctly  the 
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application. 


were  corrected  d£ 


using 


instantaneous  imi 


it,  tor  example,  t/.,a  ori 

les,  the  magnitude  of  the  minimum  value  was  found  to 
be  only  72,5  percent  of  the  maximum  value  necessary  in  the  worst  case 
for  launch ,  However,  if  the  corrections  were  applied  at  14-day  intervals, 
based  on  the  lunar- cycle,  the  minimum  value  is  only  63,5  percent  of  the 


maximum  value.  These  results  should  play  a  very  important  role 


ultimate  design  of  the  thruster  assembly  and  impulse  programming  control 
system  aboard  the  satellite. 


The  results  of  this  study  indicate  that  further  study  should  be 
applied  to  the  optimization  of  the  magnitude  of  the  corrective  impulse 
necessary  by  l)  variation  of  the  orbital  parameters  about  their 


nominal  values  ana  z;  oy 


2)  by  de termini! 


of  an  ot 


impulse 
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program  timing  for  the  thruster.  Study  could  also  be  extended  to  a 
similar  analysis  of  orbits  which  have  periods  of  revolution  which  axe 
integral  fractions  of  the  sidereal  rotational  period  of  earth  about 
its  axis. 
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»ndix  A 


Perturbation  Equations 


Oblate  Earth,  The  most  general  expression  which  is  used  to  describe 
the  gravitational  potential  of  the  oblate  earth  is  given  by  Baker  (Ref  3« 
174)  as 


,  .  !k  l+  [v«(^  ♦  s^oJl  W 

r  n*2Lk»o\  r  /  J, 

The  Legendre  polynomials  (sin  6)  and  the  harmonic  coefficients  Cnk  and 
axe  developed  and  tabulated  at  the  end  of  this  section.  The  east 
longitude  of  the  vehicle  is  given  by 


kB  "  ^  “  °go  "  %  ^  “  V 


The  harmonic  function  given  above  can  be  separated  into  four  ports  and 
Bivon  h,  *  -  ♦  *3  "  .  man  tte  tvo-toiy  un  Is 


The  reeatning  terns  are  the  zonal  harmonies* 


v  t> o 


tf  tala  6? 


where  J  are  the  zonal  harmonic  coefficients  tabulated  later?  the 
n 

tesseral  harmonics, 
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and  the  sectorial  harmonics 
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The  contributions  to  gravitational  acceleration  c*  a  body  in  earth 
space  is  given  in  terms  of  the  gradient  of  the  potential  function  as 


x  ♦  r 


-  v<Px  *  V<P2  +  +  V^4 

where  the  two~fcoiy  acceleration  vector  is  given  by 


(20) 
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where  the  two-body  posit-ion  ecapoaeats  auat  be  these  of  the  reference 

»V 

orbit.  The  perturbative  contributions  r  *  due  to  the  &ml,  tessera! 
and  sectorial  Jarccnics  axe  given  by 
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where  (sin  6)  baa  been  replaced  by  its  equivalent  , 
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where  vX„  «  Va  -  ^(tan  ^  ),  which  gives 

JCr  X 


dot  -y  x  dot 

_  .  ___  ,  _  .  «d  yj  -  0  (25) 

the  adopted  zonal  haxtuonie  coefficients  as  used  la  the  analysis 
are  given  in  Table  II  (Ref  9t0). 


table  Hi  Zonal  Harsonlc  Coefficients  * 


8 

Ja 

2 

10S2.&S  x  ID*4 

3 

-  2,536  x  lO'6 

4 

-  1,593  x  10*4 

5 

-  0*230  x  104 

6 

0.502  x  10*4 
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The  normalized  tesseral  and  sectorial  harmonic  coefficients  used  in 

viie  analysis  are  presented  in  Table  III  (Ref  9*60),  with  the  normalizing 

coefficients  N  .  obtained  from  the  relation 
nk 


[2  (2n  +  1)  (n  -  k)M 


(n  +  k)i 


(26) 


Table  III,  Normalized  Tesseral  and  Sectorial  Coefficients 


r . . 

‘V 

a* 

k 

V  (Jd°6) 

5nk  ^ 

"nS 

2 

2 

2.41290 

-1.36410 

0.64550 

3 

1 

1.96980 

0.26015 

1,08010 

3 

2 

0.89204 

-0. 63468 

0.3415? 

3 

3 

0.68630 

1.43040 

0.13944 

4 

1 

-0.52989 

-0.48765 

0.94868 

4 

2 

0.33024 

0.70633 

0,22361 

4 

3 

0.98943 

-0.15467 

0,05976 

4 

4 

-0.07969 

0.33928 

0,02113 

5 

I 

-0.05362 

-0.09791 

0.85635 

5 

2 

0* 61286 

-0.3503? 

0,16183 

5 

3 

-0.43083 

-0.08686 

0,03303 

5 

4 

-0.26693 

0. 09301 

0.00779 

5 

5 

0.12593 

-O.59910 

0, 00*46 

6 

1 

-0,09898 

0.03765 

0.78680 

6 

2 

0,05482 

-0,35175 

0.12440 

6 

3 

0.02787 

0,04483 

0,02073 

6 

4 

-0,00040 

-0,40388 

0.00379 

6 

5 

-o,ai43 

-0,52264 

0,00081 

6 

6 

0.08869 

-0*07476 

- - - 

0,00023 

The  values  of  ti»  conventional  coefficients  (C^  and  S^)  are 
determined  by  the  relations  and  . 
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The  Itf'geadre  polynomials  which  are  required  in  the  zonal  and  tesseral 
harmonic  functions  can  be  obtained  from  the  Rodrigues  formula 


(x) 


1  dn  (x2  -  l)n 
2n  a*  dx11 


(2?) 


where  Pn(x)  is  the  polynomial  of  order  n,  and  x  »  sin  6  -  -  .  The 
expressions  for  n  -  0  through  n  -  6  axe  obtained  as  given  below, 


PqOO 


(28) 


-  x  (29) 

?2(x)  ~  |  (3*2  -  1)  (30) 

P3(x)  -  |  (5*3  -  3x)  (31) 

PgU)  ■  |  (3 5x4  -  3 Ox2  +  3)  (32) 

*$(*)  *  |  (63x5“70p  +  15x)  (33) 

P6(x)  -  (23lx6  -  3l5x4  +  I05x2  -  3)  .  (3*) 


the  associated  Legendre  polynomials  are  obtained  using  the  relation 


*£(*}  -  u  -  *2) 


k/2  4k  r *(*) 


dx 


D5) 


where  P®(x)  •  £.,(*)♦  The  expression  which  result  are  gives  as  follows} 


r^(x)  »  0  ,  n  <  k 

(36 

?*(x)  -  (1  -  x2)1^2 

(37) 

pj{x)  -  3x  (1  -  **)V* 

0») 

P^U)  -  |  (5x2  -  1)  (1  -  X2)1/2 

(39) 
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rjw  ■ 

\  (7X3  -  3x)  (1  -  x2)1'2 

(40) 

I5W  - 

^  (21x4  -  l4x2  +  1)  (1  -  x2)1/2 

(41) 

t\(x)  = 

y  (33x5  -  30X3  +5x)  (1  -  x2)1/2 

(42) 

4m  - 

3  (1  -  x2) 

(43) 

4M  • 

15*  (1  -  x2} 

(44) 

$*>  - 

^  (?x2  -  1)  (1  -  x2) 

(45) 

•jW  • 

^  (3X3  -  x)  (1  -  x2) 

(46) 

4.M  • 

2~P  (33X4  -  lgx2  +  1)  (l  -  x2) 

(4?) 

I^(x)  » 

15  (1  -  x2)3/2 

(48) 

P4U)  ■* 

105x  (1  -  x2}3/^2 

(49) 

P'jcx)  - 

^  (9xL  -  1)  (1  ~  x2)3^2 

(50) 

p|(x)  - 

(u*3  -  3x)  (1  -  x2)3'/2 

<51) 

*£(*}  - 

105  (I  -  x2)2 

(52) 

»jw  - 

9i$z  (1  -  x2)2 

(53) 

»*«*>  * 
f|(x)  » 

«  (1  .  x2)2 
& 

945  (1  >  x**r/fc 

(54) 

(55) 

?|(x^  * 

103$3x  (1  -  x2}5/2 

(55) 

4<*>  * 

10395  (1  -  x2)3 

(57) 

The  caajseaents  of  the  gradient  of  the  associated  Legendre  polynomials 

k 

given  generally  by  W 

and  expressed  in  equatorial  coordinate- 

a 

are  deterainsd  fro©  the  following  re la* ionsj 
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1  t  (?) 

=.  A  (S:)  A-  (—] 

bx  \r /  ^zj  n  \r/ 

(58) 

A  pk  /s' 

by  n  \r/ 

«  A  (s\  .  JL  A  i-\ 

by  \r /  j  n  \r/ 

(59) 

Hit). 

a  At*).  *  pk/*) 
da  lr/  ^sj  n  \r/ 

(60) 

o 

where  r*- 

-  xZ  +  y2  +  s2 

(61) 

£(9 

»  xs 

V 

(62) 

A  /s\ 

by  irj 

-  .  la 

P 

(63) 

0%  (r) 

2  2 

p 

(64) 

J-L  ^  n\ 

.  -Ha ,  .  $$  (?) 

a(?)  »u-) 

[*  -0^  1  -(Jf 

(55) 

Substitution  of  Eqs  (6a)  through  (S3)  lato  &JS  (36)  through  (&G)  will 

give  the  required  coapsaesis  of  the  gradient  fxmcticas. 

THlrd  Body  .^traction,  Hie  literal  equation  oZ  aotl,©»  lor  si  snail 
Basis  such  as  aft  earth  ©stall  11;  is  the  presence  of  a  a&  la  attract  lug  body, 
the  earth,  as  veil  as  other  poriurc-isxg  bodies  is  given  tsy 


««  *\ 
••  •'  '(  I1! 

K  jX] 

(  l 

v-*  „ 

[i(xl  1 

k  *  —  ■ 

•r3  ' £  * j  ?pj 

x ' 

r  +  r  *  i?  ’  » 

hi  Is 

’  'P|y 
u 

-  / 

w.  J 

j 

y 

£ : 

m 


!x 
y 

p. 


*  radial  vector  between  the  ire  in  and  port.  :rhing  tody 
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“  radial  vector  between  the  satellite  and  j  body 


'  *5  *  4  +  4  *  4 

4 "  4s*  4*  4*  (x  ■  */* (y  -  ?/* (z  ■  92  (68} 

The  first  term  of  Eq.  (66)  again  represents  the  simple  2-body  attractive 
acceleration  while  the  second  represents  the  perturbative  acceleration 
due  to  the  third -bodies  (Ref  i:3 89). 


The  perturbative  equations  for  the  case  of  lunar  and  solar  masses 
can  be  written  as 


i  [*'%) 

rTly-y* 

pa  U  -  zj 

i  \*-s- 

~  y"  V 

r 

ps  a  -  .?-  J 


l  ®  A. 


‘e  Is]  J 
s  ■ 


The  radial  components  x  ,  y  and  are  available  directly  from  the 
ephemeris  tables  (Kaf  3  0^)f  although  the  components  and  a  must, 

bs  calculated  using  the  tabulated  right  ascension,,  declination  and  semi- 
diameter  fir  the  moon  (Ref  Is 52)*  These  calculat ions  can  be  mad®  by 


using  the  relations 


x  *»  r  cos  a 

ms  ® 

y„  “  r  sin  ct 

vju  m  ® 

a  “  r  sin  6 

m  m  ro 

r* 

rm  “  s 

m  e  « 


where  r'  is  rho  reference  distance  given  a  semi-diameter  s'  , 
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A^giin, the  components  x,  y  and  z  must  represent  the  position  in  the 
reference  orbit  as  shown  by  the  vector  a  in  Fig,  29  . 


Vehicle 


Figure  29.  Third -Body  Attraction  and  Solar  Shadow  Geometry 


Solar  floatation  Pressure ,  Consider  a  simple  case  of  a  spherically 
symmetric,  perfectly  absorptive  satellite  of  mass  m  in  earth  orbit.  The 
perturbative  acceleration  caused  by  the  radiation  pressure  of  solar 
photons  impinge nt  on  the  surface  area  A  presented  to  the  photons  is 


i  [*-v 

po  s  rr y  -  ys 

rpS  2  -  2 
r  B 


where  f  is  the  fraction  of  the  solar  surface  presented  to  the  vehicle 
s 

as  the  vehicle  passes  behind  the  earth  and  Pq  is  the  solar  constant, 

A 

Baker  has  indicated  that  for  area-to-mass  ratios  -  which  are  ss  than 

m 

4f\  - 

ID  cnr/gm  the  calculation  need  not  include  a  variable  factor  fg  ,  that 
is,  a  point  source  sun  may  be  assumed  (Ref  2s 186-197).  In  this  case, 
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f  “  1  if  the  vehicle  is  in  the  sunlit  portion  of  the  orbit  and  f  *  0 
s  s 

if  it  is  in  the  shadowed  areas.  Referring  to  Fig,  29 t  the  factor  f  can 

s 

be  obtained  by  using  the  relations 


with  f  “  1 


^A 


cos 


-1  /rT>S*' 


a  (>r  ®  f  -  o  if  f 


0  (for  point  sun) 


<  y. 


(75) 

(76) 

(77) 
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The  radial  distance  r  in  the  orbital  plane  is  given  by 

P 


1  +  e  cos  v 

where  the  true  anomaly  \)  is  given  by  the  relations 

cos  E  -  e 

cos  V  =  - - 

1  -  e  cos  E 

a  '/l~-~e^ 

sin  v  *»  — - - —  sin  E 

r 


(83) 

(84) 

(85) 


In  terms  of  the  perifocal  coordinate  system  PQW,  the  radial  and 
velocity  vectors  are  given  by  the  relations 


r  * 

A 

r  cos  v  P  + 

r  sin  v  Q 

(2) 

y  is  - 

A 

\  [  “  sin  y 

V  P 

P  +  (e  + 

COS  V  )  Q] 

(3) 

where  P  «  unit  vector 

in  the  direction  of  perigee 

A 

Q  «  unit  vector  perpendicular  to  P  (in  the  direction  of  motion) 


A  A 

The  unit  vectors  P  and  Q  can  be  expressed  in  terms  of  the  unit  vectors 

A  A  A 

I,  J,  and  K  in  XY 7>  space  with  the  use  of 

A  A 

P  »  (cosuicosft  -  sinwsinQcos  i  )  I 

+  (coswsinfi  +  sin U)  cos  il  cos  i  )  J  +  sinwsin  i  K  (86) 

A  A 

Q  “  (-oosS?  sinu  -  sinftcoswcos  i  )  I 

/  A  A 

+  (-sinwsin^  +  cosocosft  cos  i  )  J  +  ooswsin  i  K  (8?) 

Substitution  of  Eqs  (8b)  and  (8?)  into  Eqs  (2)  and  (3),  using  the  relations 

cos  u  ■  cos  (<j  +  y)  •*  cosu) cosy  -  sinu^siny  (88) 

sin  u  «*  sin  (w  +  y)  «  sinwcosy  +  coswsinJ'  (89) 
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give  the  position  and  velocity  components  in  the  XYZ  coordinate  system! 


x 

y 

z 


r  (cosficos  u  -  sin  ft  sin  u  cos  i  ) 
r  (sin£2cos  u  +  cosSlsin  u  cos  i  ) 
r  sin  u  sin  i 


(90) 

(91) 

(92) 


K 


-  sin  v  n|-—  (coswcosft  -  sinwsinft  cos  i  ) 

Jr 


'Ke 


+  (e  +  cos  v  )  -jf-—  (-cosftsintu  -  sinftcoswcos  i  ) 


(93) 


K 


-  sin  v  i|“  (cosn/sinSl  +  sine;cosQ  cos  i  ) 


K 


+  (e  +  cosy)  i|~  (~sin£2 sinu>  +  cos£2coswcos  i  ) 


K 


-  sin  v  -J~~  (sinwsinl)  +  (e  +  cosy  (coswsin  i  ) 


(94) 

(95) 


The  total  orbital  velocity  is  given  by 


(  i2  +  y2  +  i2  )V2  (96) 


Osoulated  Orbital  Elements.  Given  the  set  of  position  and  velocity 
components  at  any  point  in  an  orbit,  a  set  of  orbital  elements  can  be 
calculated.  Equation  (96)  car  be  solved  for  the  semi-major  axis 

-1 


a 


2 

,r 


K 


e 


(9?) 


The  components  of  the  orbital  momentum  per  unit  mass  h  are 


hx- 


yz  -  ay 
zx  -  xz 


Kith 


z 

h2- 


xy  -  yx 


(98) 

(99) 

(100) 


*1  *  hy  +  4 


(101) 
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The  semi-latus  rectum  p  and  eccentricity  e  can  he  determined  hy 

h2 

p  ■  r 


d-f)1/2 


(102) 

(103) 


The  eccentricity  vector  and  its  components  are  given  as 


«  -  i £(v‘-^)F 


K 


(r  *  v)  v 


ex 

1  , 

K 

X 

1 

X 

ey; 

,  o 

—  (v2 
K  K 
e 

.-a-). 

r  y 

y 

-  r 

• 

y 

W 

,z. 

z 

r  •  v  **  xx  +  yy  +  zs 


(104) 


(105) 


(106) 


The  balance  of  the  elements  can  now  he  calculated.  The  inclination 
is  given  by 


cos 


-1  /V 


) 


(10?) 


The  longitude  of  ascending  node  is 


ft  *  cos 


-1 


*h 


hfc  +  h2 
x  y 


(ioa) 


where  if  h  <  0,  ft  >  180  ,  The  argument  of  perigee  is 

X 

,  /  *h  e  +he  , 
w  -  cos  -1  ( 

W*4  +  hv 


(109) 


where  if  <  0,  w  >  180  .  Finally',  the  true  anomaly  is  given  by 


.1  /  V  +  v  *  v 

V  m  cos  . . . . . . 

\  er 


) 


(no) 


where  if  (r  •  v)  <  0,  v  >  180  , 
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Impulsive  Correction  of  Perished  Orbit,  Very  useful  relations 
have  been  developed  and  tabulated  by  Ehricke  (Ref  ?i35^)  and  axe  applied 
in  the  present  section.  Attention  must  be  directed,  however,  to  the 
basis  of  their  usefulness.  The  impulses  which  are  calculated  using  these 
relations  are  only  accurate  if  they  axe  much  smaller  than  the  velocity 
components  of  the  vehicle.  It  is  further  assumed  that  the  impulses  are 
made  instantaneously.  The  geometry  of  the  orbital  correction  is  shown  in 
S.  30  . 


The  inclination  and  the  longitude  of  the  ascending  node  can  be 
corrected  simultaneously  by  making  a  plane  change  with  an  orthogonal 
impulse  whose  magnitude  is 


Ai  Vj  AO  v,  sin  ij 

■A  V  " 

(12) 

"  cos  (vL  -  t/A)  Sin(^-PQ) 

* 

a  a 

(l U) 

at 

\  "  %  *  tan  XI  sin  V 

where  is  the  true  anomaly  at  the  ascending  node  and  b  represents 
the  point  at  which  the  correction  is  made,  The  axgrasant  of  perigee  is 
corrected  at  either  the  apogee  cr  perigee  with  a  normal  impulse  of 


Av, 


*A  *A  VA 


ft 


^A  +  rA  C“  {\  -  " 0 > 


(«; 


wtiere  i>Q  is  the  true  anomaly  at  perigee  and  the  subscript  k  represents 
values  at  apogee.  This  normal  impulse  will  result  in  a  negligible  change 
in  tlw  eccentricity 


A  e 


r  dv  sin 
A  „  a,..- . A. 

SA  VA 


since,  v '  ^  180  and  ein  ^  ~  0  * 


(112) 
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Figure  30.  Geoasiry  of  Oreat-al  faraxaster  Correction 


The  tangential  iapuleo  aadt'  at  apogee  la  determined  by  tho  value  of 
chains  required  to  correct  tho  perigee  distance  r^  »  and  is  related  to 
the  change  required  in  the  seal-major  axis  as  follows 


A*A  * 


r  -  r 
PQ . *P 


(113) 


Tho  required  ixspU.ce  is  found  to  be 


ivu  - 


aa.< K- 

2 

A  aA  VA 


(W) 


which  results  in  significant  changes  in  tl»  eccentricity  and  Bean  notion 
as  given  by 
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Ae,  .  *»  2  (e.  +  cos  v, ) 

■  a  A  A 


Av 


ta 


AM.  «  - 
A 


3  VA  Avta 
KeaA 


(114) 

(115) 


la  addition,  a  negligible  change  results  in  the  argument  of  perigee  as 
compared  to'  that  due  to  the. normal  impulse.  This  change  is  given  by 


Ao,  .  2  C'A  -  vo)  AVU 

A 


®A  VA 


(116) 


Upon  arrival  at  perigee,  a  second  tangential  impulse  is  required  to 
fully  restore  the  values  of  the  semi-major  axis  and  eccentricity  to  the 
nominal  values.  The  increments  of  change  required  are 


A°P  •  ao  -  ("r  +  4»a) 

A.p  .  v  -  (&p  +  iaA) 


(11?) 

(118) 


The  velocity  at  perigee,  prior  to  the  impulse,  is 


'  1  ■ 

1 

Ke 

,rP0 

(ap  +  Aa^J 

1/2 


• (119) 


and  the  eccentricity  and  semi-major  axis  valws  arc 


ei  *  ftp  *  A*k 


*p  *  ap  + 


(120) 

(121) 


The  tangential  impulse  may  now  be  calculated  using  the  relation 


Av 


Ae  v*  Aa  K 

V  P  „  p  e 

tp  "  2  (o^  +  cos  ITT  2  ^ 


(li>) 
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and  will  result  in  changes  of  the  mean  motion  and  argument  of  perigee 


n 


AO) 

P 


3v 


2  sin  ( V  -  v)  Av. 

P  o  tp 


e  ’  v* 
P  P 


(122) 

(123) 


The  change  of  mean  motion  is  of  no  concern  to  the  analysis  since 
integration  of  the  succeeding  orbit  will  proceed  from  this  point.  The 
change  in  argument  of  perigee  will  be  vanishingly  small  since  the 
nominal  value  will  have  been  achieved  by  the  normal  impulse  at  apogee. 


The  values  of  the  orbital  elements  (essentially  the  nominal  values) 
which  are  obtained  after  the  corrective  impulses  are  given  by 


a 

p 

ap 

+ 

Aa.  +  Aa„ 

A  P 

(12*0 

a 

m 

e 

P 

+ 

Ae  +  Aeyv 

A  p 

(125) 

i 

sa 

i 

P 

Ai 

(12^) 

(i) 

- 

u 

P 

A4> 

(127) 

n 

• 

n 

TJ 

+ 

Ail 

(128) 

The  time  of  arrival  at  perigee  after  applying  the  impulses  is 
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Appendix  C 


Numerical  Techniques 

Encke's  Method,  One  of  the  more  refined  techniques  available  for 
computation  of  orbital  perturbations  was  developed  by  J,  F.  Encke,  The 
principal  advantage  of  using  this  technique  is  that  only  the  perturbative 
accelerations  need  be  integrated.  This  leads  to  greater  accuracy  while 
using  larger  integration  steps  as  compared. to  the  Cowell  method,  in 
which  the  sum  of  all  accelerations  axe  integrated  together  (Ref  8» 225, 
6i469). 


The  equations  of  motion  for  the  perturbed  and  reference  orbit  shown 
in  Fig,  31  ore 


r* 

r 


** 

P 


K  _  * 

•4r  +  r 

Ke  - 

w  0 
r 


P3 


(130) 

(m) 


respectively,  whore  r  is  the  sir®  of  all  the  perturbative 
accelerations  diecuo&sd  in  Appendix  A, 


Figure  31,  Ceoetry  of  Perturbed  Orbit 
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Considering  only  the  x-component  for  purposes  of  illustration, 
noting  that  6r  *=  r  -  p  ,  and  taking  the  difference  between  the 
perturbed  and  reference  orbital  equations  (130)  and  (131) ,  we  liave 


#  *»  •» 
6  x  “  x  -  p 


K 


K. 


A_ 

l 

x  -  6x 


x 

*1 

X' 


3 

J 


.V 

+  X 


*v 

+  X 


1  -  £  x  -  6* 


+  X 


(132) 


Since  we  also  have  the  relations 


r2  «  x2  +  /  +  zZ  «»  (P  +  fix}2  +  (p  +  6y)2  +  (p  *  5s)2 

x  y  s 

“  PZ  *  2Px6x  +  2P  6y  +  +  6x2  +  6y2  +  6s2  (133) 


we  can  write 

£  .  1 


1 


(pv  +  §&x)5x  +  (p  ♦  t$y)6y  +  (P  * 


V  2 


¥  * 


P* 


(13*0 


Defining  the  quantity  in  the  bracket©  as  the  dlss&nsionlesn  fa&lor  q, 


P* 


1  ♦  2q 


and 


p- 


-3/2 

*  (1  +  2q)  U35) 


which  can  now  be  used  to  rewrite  bq  (132)  as 


6x  “ 


K  J 
e  1 


3 

P*'  t  L 


-3/21 

“  (1  +  2^5  J*  - 


5xj  ♦  x  (136) 


Expanding  the  torn  (1  +  Sj.)  '■' 2  in  a  bineaiLa.1  series  gives  tbs 
re  aidts 
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fq 


i 

3q 


-3/2 

■  (1  +  2q) 

3*5  2  3v5*?  .3  3*5;7*_9  4 

-  *r±  q  +  vT  q  -  4T  4 


25  "  31 


3.5*?*9*X1.5 
*  — 5i  q 


(13?) 


6x 

K 

“ 

X 

■fix’ 

X 

.tv 

6y 

e 

w  — — 

„3 

fq 

y 

- 

6y 

+ 

y 

•v 

Az 

PJ 

.  a 

6a 

a 

(13s) 


The  integration  of  tic,  for  example,  v>  determine  6x  requires 
the  knowledge  of  6x  and  x  ,  aid  In  therefore  an  iterative  process. 
Tbs  first  application  of  tht  w  orical  integration  of  %  (133)  can  be 


based  on.  6x  **  0  and  x  -  «^i c‘-»  results  in 


to 


K  P 
"ex 


♦V 

*  x 


(139) 


«vi  tho  flni  .bk*K»««  «'  »*•  Th»  “cond  lt»r*Uo"  tosl'w 
IMofpaUon  of  2%  (138)  «W>  *  *  »x  *  **•  t!w  lt*“U<*  U 

o«itUw»4  until  luoisalfieont  chsagoa  »»  «“»*  lfe<l  s*' 

tt*  mittU  In  this  theeis  lndlcUd  th»t  sin  ItSfsUoho  «« 

.uKUlsnl  for  this  par-.,  *>>»«  *  «»»*««  *“*  o£  •»  {U?)' 
iooortfc*  to  Ehrleta  (*rf  6<h?5i.  tiw  «*»»««»  ittaUon* 

Involve  roploccfiont  of  Uth  «*“»  rf  1  ’**'*  *°*W 

slso  topu  the  m->SW««l  ‘®  r.«0«l»U  the  values  of  X  »»*»  «* 

lstor  nlw  Of  *  ntlwr  than  2»  (Kbit  »..»  Ihtogsetod  using  «*• 
nocuous,  ifter  osch  of  M*>  »  sot  cf  oocuUUU  ubtul  •!— *•  m. 
calculated.  »  »  Olosonts.  cor«oton41«  to  U»  S»>»t  of  soon  Motion, 

uoro  then  used  to  ««.»«  «*  *«  Section  «•*  «»  »*««•  -hit. 
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Method  of  Numerical  Integration,  The  step  integration  of  the 

pertur tat ions  was  carried  cut  by  application  of  the  Gauss -Encke  met!  od 

of  mechanical  quadrature  (Ref  6:486),  This  technique  is  also  called  the 

2 

Gauss-Jackson  or  Sum-Squared  ( )  method  (Ref  4s4l8),  This  numerical 
technique  requires  the  construct!  jf  a  table  of  data  based  on  equally 

p 

spaced  perturbative  functions  f^  -  w  6x  where  w  is  the  time  interval 
between  the  data  points,  and  6x  is  cue  perturbative  acceleration. 

Since  the  dimension  of  the  function  is  length,  all  entries  in  the 
table  will  be  in  units  of  length.  The  general  construction  of  the  table 
is  illustrated  in  Table  IV  on  the  following  page.  Note  that  the  entries 
are  tabulated  as  a  function  of  the  normalized  time  unit  t  »  t,/w  , 
where  t^  is  the  time  measured  from  some  reference  point  (e,g.,  the 
osculation  point). 


Follovfing  the  tabulation  of  the  perturbative  functions  f^,  the 

first  step  in  the  generation  of  the  data,  table  involves  the  calculation 

of  the  differences  6™  using  the  following  relations; 

J 


*1 

V*  " 

f 

(140) 

- 

6i 

.n-1 

h+i  ■ 

An-1 

h-i 

(141) 

,  n— 1  „ 
6i+i 

.n-1 

6i+i  • 

5„-i 

(142) 

where  i  a  ,-3,-2, -1,0, 1,2,3* . ♦ 

,  ,  and 

n  "  2,3 1 ^i5»  •  * * • 

Assuming  that  the  higher  order  differences  approach  a  constant  which  is 
very  small  compared  to  the  function  f^,  the  differences  in  the  upper 
a..d  lor?  right  hand  corners  of  the  table  can  be  approximated  by  using 
Eqs  (l40'i  through  (142),  The  present  analysis  considered  as  constant. 
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Table  IV,  Gauss-Sncke  Hnaerical  Integration  Tabular  Construction  ■ 
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6x 
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6x_3 

wox^ 

Y2 

-2 

6x-2 

*6i-2 

I2 

-2 

-1 
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«6x^ 

E-1 

0 

6x0 
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1 

6xl 

"6*1 

K 

2 

6X2 

•.«>x2 

3 

6X3 

w6x^ 

0 

e: 

3 

4 

6*4 

“5x4 

E* 

4 

I  V 


^  f 


'2  c2 
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I  -i 


**-&  ,2 


e; 

-v2 


"2  f. 


f2  f. 


* 


r/2  f„ 


6*  '  *  64 
0»1  0_1  - 

62  *  ‘-fe 

0  1  .3  ■  0  1  ,5 

W  V 


%  2  %  '  4  \ 

1  6J  I  >2  .«£  ^ 

I  A,.  2 


1  2 

The  calculations  of  the  suras  E  and  E  0X6  not  as  straight¬ 
forward  as  those  of  the  differences.  Application  of  central  difference 
integration  formulas  gives  the  relations 


‘<r 


9  1  1  1  11 

26xdr  «  (E,»)  -  —  (6^)  +  — 


(si>  -  5  (6i> +  ts  <si> 

191  5  2497  7 

60480  +  3628800 


(143) 
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6xi  =  Jj* !?'  6x  d7^ 


Ei  +  12  fi  “  240  6i  *  60480  6i 


3*  ,4 


289  6« 

31? 

(144) 

36288OO  i 

22809600 

,  1. 

& 

/  1 

1  \ 

( T  \ 

a  *— 

(l  .  <  “ 

£,.a  ) 

(145) 

'  % 

O 

V  *4 

i+T/ 

<«*>, 

n 

« 

2 

/  - 

iffy 

(146) 

.  Aside  froia  the  perturbative  functions  and  the  difference  table 

which. has  been  compiled,  the  values  of  ox  and  w  6x  at  r  «  0  are 

o  o 

1  •  }  .  2 

also  known  and  are  used  to  calculate  £  ■  , %  ^  and  with  S^s  (143) 

.  2  ”2*  0 

and  (144),  Note,  .however,  that  since 


r*  -  r3-  -  •*„  :  m 

rj  +  “  2  K  -  fo  (1W> 

*  *  • 

we  can  rewrite  Eqs  (143),  ( 144)  and  (14?)  to  obtain  the  relations 


1 

y  « 

H 

1  It  11  0 

0  20  12  VV  720  v  o' 

191  ,„«>x  2497  ,.7 

+ - (64  -  - - —  (6  J  +  *♦••• 

60480  0  3628800  0 

(149) 

£  1  “ 
4 

s1  -  f 
i 

(150) 

„2 

£  * 

0 

1  1  2  31  4 

6xo  "  12  fo  +  240  6o  “  60480  6o 

.  289  .6  317  .8  . 

t  0 

3628800  °  22809600  ° 

(«i) 
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Once  these  three  key  values  are  determined ..  the  balance  of  the  £  and 
2 

Y,  columns  of  the  table  can  oe  filled  in  using  the  following  equations 


v 


i+|- 

2 


zl .  +  fA 

i-if  i 

(152) 

E2  +  E1 
i-l  i-f 

(153) 

Finally,  the  values  of  ox^  and  w  6x^  can  he  calculated  using 
Eqs  (143)  and  (144), 


Interpolation  Equation,  The  calculation  of  the  perturbations  due 
to  solar  and  lunar  attraction  requires  interpolation  of  the  ephemeris 
data  which  is  provided  at  daily  intervals.  Determination  of  the  times 
of  perigee  and  apogee  passage  is  also  necessary  during  the  course  of  the 
analysis  and  requires  the  interpolation  of  the  position  and  velocity 
components,  A  useful  equation  is  the  fourth-order  central -difference 
interpolation  polynomial  (Ref  10 i 216)  given  by 


V  » 


(t ) (t-w) (t-2w) (t+w ) 

'  £ 

24  w 


(t)(t-w)(t-2w)(t+2w) 

y  +  - - - -  y 

~2  -6i* 


(t-w)(t-2w)(t+w  )(t+2w) 

4 

4v? 


Y 


(t)(t~2w)(t+w)(t+2w) 

-6w4 


'+1 


( t ) ( t-w ) ( t+w ) ( t+2w ) 
24w4 


J+2 


(154) 


where  w  is  the  tire  interval  between  the  equally  spaced  data  points 
(y_2  through  y+2)  and  -2w  <  t  <  +2w. 
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Determination  of  Extrema,  The  determination  of  the  times  of  perigee 
and  apogee  passage  involves  the  location  of  minimum  or  maximum  radial 
distances  respectively  among  five  equally  spaced  points,  A  fourth-order 
polynomial  (Ref  10:223)  can  be  fitted  to  these  five  points  and  is 
generally  given  by 


y  “  V  +  &y  +  v  +  + 


(155) 


where  the  coefficients  can  be  determined  from  the  relations 


a,, 


a,,  *> 

3 


a. 


“1 


a. 


y-2  -  4y-l  +  6y0  "  %+l  + 


+2 


24w 

■y_2  +  2y_1  “  2y+1  +  y. 


+2 


1Z\P 

-y_2  +  I6y-1  -  30y0  +  I6y+1  -  y 


+2 


24s/ 


y_z  -  Sy^  +  8y+1  -  y. 


+2 


12  w 


(156) 

(157) 

(158) 

(159) 

(160) 


The  value  of  t  for  the  maximum  or  minimum  of  the  function  y  con  be 
found  by  the  solution  of  the  expression 

fj:  "  0  »  4a^t3  +  3a^t2  +  2agt  +  a ^  (l6l) 

by  using  Newton's  method  as  described  in  Appendix  B,  whore  the  function 
f  is  given  by  Eq,  (l6l)  and 

f»  «  I2a4t2  +  6a3t  +  2a2  (162) 
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Appendix  D 


Computer  Model 

General  Description,  The  computer  model  written  for  the  analysis  of 
the  problem  described  in  this  thesis  consisted  of  the  executive  program 
PROGRAM  ENCKE  and  appropriate  subroutines  which  were  used  to  make  repeated 
calculations  during  specific  phases.  This  section  is  devoted  to  a  brief 
explanation  and  description  of  each  segment  of  tne  model, 

1)  PROGRAM  ENCKE  is  used  to  provide  '  the  source  of  input  of 
initial  parameters  and  ephemerls  data  and  output  of  desired  data  following 
the  various  phases  of  analysis.  It  also  provides  the  framework  for 
relating  the  calculations  performed  by  the  subroutines.  The  listing 
presented  in  the  following  section  specifically  represents  the  case  which 
includes  daily  correction  of  the  perturbed  orbit,  A  flow  diagram  which 
proa-ides  for  other  cases  is  shown  in  Fig,  32. 

2)  SUBROUTINE  POSIT  calculates  the  position  and  velocity 
components  of  the  reference  orbit  following  the  generation  of  a  set  of 
eccentric  anomaly  values  by  solution  of  Kepler's  equation,  Since  the  values 
are  equally  spaced  at  small  time  intervale ,  successive  trial  values  of  E^ 
are  assumed  to  be  those  of  the  point  just  determined.  Values  for  all 
orbital  points  are  generated  at  a  single  call  from  PROGRAM  ENCKE, 

3)  SUBROUTINE  POLYLEG  calculates  the  values  of  the  Legendre 
polynomials  and  their  gradients  for  each  orbital  point  as  called  for 
individually  by  PROGRAM  ENCKE. 

4)  SUBROUTINE  ZONAL  calculates,  the  value  of  tho  aonal  harmonic 
contribution  to  the  oblate  earth  portui. Native  acceleration  for  each  orbital 
point  as  called  by  PROGRAM  ENCKE. 


5)  SUBROUTINE  TESSER  provides  the  tesseral  contribution  of  the 
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Figure  32.  Flow  Diagram  of  Computer  Kodel 
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oblate  earth  perturbative  acceleration. 

•  6)  SUBROUTINE  SECTOR  provides  the  sectorial  harmonic 

contribution  of  the  oblate  earth  perturbative  acceleration, 

7)  SUBROUTINE  LUNAR  computes  the  third-body  attractive 
perturbation  due  to  the  lunar  macs  for  each  orbital  point  as  called  by 
PROGRAM  ENCKE, 

8)  SUBROUTINE  SOLAR  computes  the  third-body  attractive 
perturbation  due  to  the  solar  mass  for  each  orbital  point  as  called, 

9)  SUBROUTINE  RADIATE  calculates  the  perturbation  due  to 
solar  radiation  pressure  for  each  orbital  point  as  called  by  PROGRAM  ENCKE. 

10 )  SUBROUTINE  TABLE  calculates  the  components  of  6r  and  vr6r, 

given  the  component  values  of  the  acceleration  functions  w^6r.  The 
integration  technique  is  the  £  -  method  discussed  in  Section  2  of 

Appcnuix  C.  The  x,  yf  and  z  components  are  computed  separately  in 
successive  calls  from  PRO-CRAM  ENCKE. 

11 )  SUBROUTINE  IMPROVE  applies  Snoko ' s  technique  to  improve 
tha  orbital  integration  initially  provided  by  TABLE.  The  orbit  will  be 
improved  using  the  number  of  iterations  given  by  the  in+egor  NN,  TABLE 

successively  called  to  integrate  the  components  of  the  improved 
perturbative  acceleration  functions, 

12)  SUBROUTINE  05CULAT  calculates  the  orbital  elements  of  the 
perturbed  orbit  given  the  position  and  velocity  components  at  the  point 
of  interest  defined  by  the  inteser  N03C,  The  basic  squat  ions  used  are 
disousseed  in  Section  2  of  Appendix  B, 

13)  SUBROUTINE  I TRICES  calculates  the  position  and  velocity 
components  and  the  time  of  arrival  at  perigee  (or  apogee)  by  interpolation 
between  equally  spaced  data  points  surrounding  the  approximate  orbital 
point  defined  by  the  integers  NI’ER  or  NAPO.  The  actual  time  of  arrival 
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is  calculated  by  using  the  fourth -order  technique  described  in  Sections 
3  and  4  of  Appendix  C, 

14)  SUBROUTINE  PELT  V  calculates  the  values  of  the  instant¬ 
aneous  impulses  required  to  correct  the  perturbed  orbital  parameters, 
according  to  Section  3  of  Appendix  B, 

Computer  Program  and  Subroutine  Listings,  The  following  section 
provides  a  listing  of  the  executive  program  compiled  to  provide  daily 
correction  of  the  orbital  parameters.  The  general  subroutines  discussed 
in  the  previous  section  are  also  listed,  A  glossary  of  the  major 
variables  used  in  the  programs  is  provided  in  Section  3  of  this  appendix. 
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PROGRAM  ENCKE  (INPUT  ,  CUTFU  T) 

COMMON/BLK A/3HG( 3, 60 )  ,  ELEH  <io)  ,PARA(  10)  ,E(6Q) 
C0MM0N/8LKB/P  (3,  5)  ,QP<6,5)  ,DQ(  8  ,5)  ,0R(8  ,  5)  ,RR(6) 
COMMON/ OLK  0/  QJ  (5)  ,  C<  7 , 5)  ,S  (7,5  )  ,FAC(  7, 5  > 
COMMON/BLKO/PERT  (7  ,3,60) 

COMMON/BLKE/CIFF  <10,60)  ,TA0(4»  60) ,  H(  3 , 5 )  , T ABS  ( 2 ,7 , 3  ) 
COMMON/ 8LKF/KRS(  3,60)  ,RRM(  3,60  ) 

COMMON/ B L KG/ RE F(  8,60)  ,  EREF  (60)  ,  ELREF  (10) 

OIMENSION  RM(3,5),RS(3,5) 

81  FORMAT  (IX  ,5£18. 1Q,2F12.8,  F12,6) 

94  FORMAT  (F15.0) 

95  FORMAT  (5I4,2F15.0) 

97  FORMAT  <6F 12  *  0) 

102  FORMAT(3X,10Ei2.5) 

REAO  94,  (BJ(I),I=i,5> 

READ  94,  (  (C  (I,J)  ,  J=l,  5)  ,1=1,7) 

REAO  94,  <  (S(I,J>  ,  J=l,5)  ,1=1,7) 

REAO  54,  (  (H(I,J)  ,J=1,5)  ,1=1,2) 

REAO  94,  (PAfiA(I)  ,  1=1,10) 

READ  94,  (ELEM(I),  1  =  1,10) 

REAO  95,  M0AY,N3AY ,NN,N,NOSC,W , THETA  0 

REAO  97,((  (R I'd,  J)  ,1  =  1,3),  (RS(  I  ,J)  ,1  =  1 ,3) >  ,J=2 ,5> 

00  14  1=1,8 
00  14  J=  i, 5 
14  P(I,J)=0. 

ALPMAlsO. 

ALPHA  2=0* 

u«o. 

USAVE=0. 

RH0(l,i)=2. 64429 
RHO  <5, 1)«3, 470844 
DO  46  1*1,10 
EREF  (I)*G. 

48  SLREF(I)  =£LEM(I) 

ELREF  (3)  =ELREF  (3  )  /PARA  (3) 

ELREF (4) =LLR£F (4) / PARA (3) 

ELREF  (5)=ELR£F(5)/PAfiA(3) 

00  400  ITOSC=MDAY, noay 
AtOMAXs-360. 

00  $0  J=  i,  4 
OG  50  1=1,3 
RH(I,  J)*RM  (X  ,  JU> 

50  RS(I, J)aRS(I,J*l) 

REAO  97,  (RM(I,5),  1=1*3)  ,(RS(I  ,5), 1=1,3) 

IF  (RM  (1 ,5  )  •  GT  #RM  ( 1,  «♦ )  )  GO  TO  51 
RM  (1 , 4 )  =RM  (1 ,4)-24. 

RM(1,3)=RHU  ,  3)-24. 

RM(l,2)=RM(i,2)-24  . 

RM(l,l)=RM(l,l)-24« 

51  NPER=-1 
NAPO=-l 
MAPO=-l 

E (1 ) =~  *32 
E(6>  =  0. 

PARA ( 3) =0* 
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U=(.TT0SC  i>-L 
00  390  I  :CT  =  i,5 
DATE-  SAVE*  ( T  ECT-1)  *0*2 
IF  ( IS.CT. GT  .  .  )  GO  TO  53 

P,  INT  8;  <ELEM(T>  *T=1,5),  RHO(  1,1)  ,RH0(5,i)  ,OATE 

53  EL  M(3>=  .EH  (3 )/ PA  R  .  v3> 
cLE  1(4)  =ELEM  <4)/PAPAU> 

ELEM5J  'EL EM  (5) /PA  PA  (3) 

CALL  PObIT  (  IT03C,  ISECT,N,  W,N0S9> 

IF  (  I^OSC.  E1"1 « 1»A N1J  »I SECTi. EQ»i>  LREF4  s  U  *RHO  <i»6» 

00  20  L  1  =  1, N 
RR(1)=RH0(1, I) 

DO  54  K=  2 »  6 

54  RR(K)=RR{K-l)*RR(i) 

CALL  POLYLEG  (I) 

T=(I“6«)*W+(ISECT-i)  *  4.3-U*24. 

IF  (T.EQ.G.)  GO  TO  120 

IF  (T.EQ.24.)  GO  TO  125 
TT=T/24. 

AAA=TT*  (TT*TT-1.  >♦  (TT*TT-4  «)/24. 

AA=AAA/ ( T T ♦ 2  •  I 
80=-<AAA*4.)  /(TT-t-1  ,) 

CC=AAA*6,/TT 

D0=®  (AAA*4.)/(TT~1.) 

EE= AAA/ { TT-2  •) 

GO  TO  130 

120  00=0. 

CC=1. 

GO  TO  127 
125  CC=0< 

00=1. 

127  AA=Q. 

BfK  0. 

£E=  0. 

130  00  160  J  =  i  ,3 

RRM(J,I>  =  AA*RMU»i)  +  38*RM(  J,2)  ♦CC*RMCJf  3>  ♦00*RM(J,4J 
1  +EC*RM(J,5) 

RRS<J»I)  =AA*RSU,  1>*G9*RS(  J,2)  +  CC*RSU,  3) ♦DO*RS < J» 4) 
1  f E£*RS ( J  ,5) 

RRS(J»!)«RRr  <J,I»*PARA  (6) 

160  CONTINUE 

RPH (1,X)=15«  *RRM  (1*1 )  /PARA  (3) 

RRM  (2  *  I )  =RRM  (2  ,1 ) /PARA  13  > 

RP.M  (*  ,  I )  =PRM  (3  ,n  *  PAR  A  (5 ) 

XH=RRH  ( 3,1  )*CCS(  RRM(  1 ,  I )  ) 

YH=RRM(3,I)*SIN<RRH(i,I)) 

RRM(  3, 1)=RRM(3,:)*SIN<RRH<2fI)  ) 

RRM(1, I) =XM 
RRM (2,1) =YH 
T=T+(USAVE+U)^24, 

AT=ATAN2(RH0  (3,1)  ,RH0(2»I)  ) 

ALPHA  =  AT  +ALPHA1 

AL=AIPHA*6  6.«C.0  0*37  5  2695*  T-THETAO 
AL0=AL*PARA(3) 

IF  4AL0.LT #-  100. )  ALPHA1=ALPHA1*2*PARA(2) 
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IF  ( AL0,LT.-18Q.>  ALO=ALO*36Q. 

PERT <  7 , 1, 1) =T 
CALL  ZONAL  (I) 

CALL  TESTER  (AL,  I) 

CALL  SECTOR  < AL, I ) 

CALL  LUNAR  (I> 

CALL  SOLAR  (I) 

CALL  RAQIAT  (I) 

200  CONTINUE 

DO  242  K=1,N 
00  241  J=l,3 
00  240  1=2,6 

240  PERT  (1,J*K)=PERT  (1  ,J,K)4P£RT(I  ,J,K) 

24^1  PERTC1, J,K)  =  PERT(1,J,K)*H*H 
242  CONTINUE 

00  245  K=1 , N 
00  245  J=1 $  8 
245  REF(J,K)=RHO  (J,KI 
00  350  I=i,l 
00  320  K  =  1,N 
00  320  J J=1 s  3 

320  PERT ( 7 , J J»  K) sPERTil, JJ,K) 

00  300  J=1 , 3 
TABS< 1  »I » J)  =  0 * 

TA6S( 2, 1 » J)  =  t  * 

CALL  TAOLE  <  1 ,  J,  N) 

00  310  K=1,N 

310  PERTH,  J,K)  =  TABU,KI 
300  CONTINUE 

CALL  IMPROVE  <X,N,M,NN> 

350  CONTINUE 

00  349  1=1, N 

T=( J~€. )  *WM  ISECT-  1)  *4,  8-U*24, 

TaTHUSAVEtU)*24, 

349  PERT ( ? , 1 »Z ) =T 
00  360  1=1,1 
00  360  J=l,3 
XXwJ+5 

TAOSU,I,J)  =  C. 

T6BS(2,I  ,J> =G. 

CALL  TABLE  { I  *  J»  NT 
GO  370  K=i ,N 

RN0(U+1,KJ  sRHCUH  ,K )  ♦  TA  B{  i  »KI 
RHOUI»K)  *RHC<II  ,K)  +  TA9C2,  K)?H 
370  CONTINUE 
360  CONTINUE 

00  380  K=1,N 

RHO(1»K)=SORT  CRMOC2,<)*RMO  (2,0  *RHO(  3,K)  *RNO  <  3 ,  K)  ♦ 

1  RH0(4,K)*RM0(4,K> ) 

RHO<$,K)=SCRTlRHOt6,K)*RK)  <6,K  )  *RHQt  7 ,  K  )  ♦RHO  f  7  ,  K)  +  * 

1  RHO  ( £  ,  K I  *RHO  (8  ,  K  H 

IF  CNAPO.GT.  3*0R.K.Lr.6,0R»K.6T.?4>  GO  TO  37l 
IF  (RHO  1 1 , K I  »GG,  RNC  <  1 ,  K- 1)  ,  A»(0  .RHO  C  ,  K)  *  CT  .RHO  C  1 ,  K*  i ) ) 
1  GO  TO  369 

GO  TO  371 
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369  NfiPO=K 

371  IF  (NPER.GT.  Q.OR.K.LT,  20.OR.K.  GT.57)  GO  TO  380 

IF  I  RHO  ( i  j  K)  oL£«  RH  0(  1 ,  K-iJ  .  ANO  .  RHO  (1 , K)  .  IT  .RHO  ( 1 , K+ 1)  ) 
1  GO  TO  372 
GO  TO  380 

372  NPER=K 
380  CONTINUE 

IF  (MAPO.GT.  C.OT.NAPO.LT.O)  GO  TO  33  3 
CALL  PERIGEE <NAPO, OUl) 

DU1=  (PERT  (7,  l,NAPO)  +  OUl*W)  /24.-USAVE 
CALL  GSCULAi  (1) 

PRINT  81,  (ELEHdl  ,1  =  1,5),  RHO(l,l)  ,RNO(  5,  I )  »DUl 
SREf'(21>atLE**(l) 

EREF  (22)  =£L£n  (2) 

EREF  <23>=ELEM(3>  /PARA  {  3) 

£R£F{24)=ELEM(£i)  /P  ARA ( 3) 

EREF  (  25) -E(  EH  (5)  /PARA  (3) 

EREF(t9>=RH0U,l) 

EREF  (2G)=RFO  <5,i) 

MAPO=l 

363  IF  (NP&R.LT .  0)  GO  TO  38?. 

CALL  PERIGEE  <NP£R  jU) 

U*(PERT17,l,NP£R) +  U*W>  /24. 

CALL  OrCULAT  (1) 

PRINT  61,  (El£H(  X)  ,1-1,5)  ,  RHO  Cl  ,1)  ,RHG(S»U  »U 
!PER=-1 

£L  E  tM  { 3 1  -  E  L  E  ‘K  2 )  /  P  A  F  A  (  ? ) 

ELEH<4)aeLE«(4»/*»AOA(3) 

ELEH(*5)  sEL£H(G)/PARAU) 

CALL  OElTAV 

PRINT  102,  (cREF(I), 1=1,10) 

DU2*<U»U$AV£)-QG1 

OU2®  1 «  /t  1  *  /0L2 ♦EREF  ( 1 1>  /F4RA  (2  )  ) 

0®USA VE  *DU i* CL2 
USAV'E  =  U 

tr  i ISECT *  EG .  5)  GO  TO  390 
385  CALL  CSCULAT  (NOSC) 

IF  <XSECT.Ml.lt  GO  TO  390 
£REF<31>«ElEMl) 

EREF  I  32)  *ELEP  (2) 

EREF  <  33)  sF*  EMU)  /PARA  (3) 

390  GCNTIM/E 
400  CONTINUE 
STOP 
ENO 
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SUBROUTING  PC5JX  (  ITOSC,  ISECT,  N  ,H,  NQSC) 

COMMON/  3LK  A  /«HO(  3,  GO)  ,  ELEM  (10)  ,  PAR  Ai  1^7  ,  £  ( 6Q> 

.IF  tITOSC.GT  .1)  GO  TO  5 
IF  (1SECT.GT.1)  GO  TO  5 
Efi)=-,3E 
£163=0. 

GO  TO  10 

5  £  <63  =£  <  NQSC) 

£ (1) =  £<  NQSC-5) 

10  00  20  1=1, N 

T=FARA(3l  +  a-e»)  * H 
H»0 

IP  (I.EQ.O.ANC.IS£CT<r.C*l.  And*  ITOSC.EQ.  1)  GOTO  13 
16  F®E<I)-T*PARA<A)-£UH{2>  ’5IN(£  (IH 

FF=l,-£UM23*CO$<  E<  I)) 

H=H+ 1 

€<I)*EU)-F/fF 

IF  (H.GE  s  1  Ofi  )  GO  TO  IS 

IF  (AOS  (F)  »L  £♦  Q.  00  00  0  uuCOO  Ci)  GO  TO  13 

GO  TO  1 6 

13  IF  (I.fn.N)  GC  TO  2& 

I F  ( I  ♦  £0«  ?»  AHC  ♦  I5£C  T*  £Q*  1 .  A  NO.I  T  Q-$C  «FQ»  1 )  GOTO  20 
EU  +  1)  at  ill 

20  CONTINUE 

21  A*COS(Ct.£Ml$U 
AA«SIK(£ienf*l) 

O^COS  <£*.£*<  3  3) 

S0s?I  N(CLFKU)>- 

V3=  SO  R  T  (  PA  R  A  (1  U  ct  EH  (  6  3  3 
CO  36  I*1,H 

0JMCOS  IE  U  I»-U£M1  7)  )/?  2)  *C0S<£  (Stl! 

RHO  '  *  ,IJ*Ei£J'{t3  /(  !.«£££*<  2i*DJ 

00®  <  E 1 1«  U  )  •  SO5*  ( l  Et-tM  (2  )  *  Eti.  E  H  ( 2  )  3.  *SI-N  i  £  U  >  31  7  AMO  (1.0 
CCC* At  AN2 ( C 0  » C ) ♦ EL  f  M  ( U 
C=COS <CCC) 

CC«$lK(lCC> 

RfcO  ( 2 ,  U  *  C *•  OM * C C* 3  J  •  F -*-0  C i,  n 
RHO ( 3  » I  >  s  (  A  A  *C *A  *CC*  U  >  *®«0  U»I> 

RNOUf  I)w  <CC*u€*>  *RK>(  I,!* 

*ri=-v3*co 
v?s(£i£H(2no)**n 
CaCOS <  £L£H  !U) 

CC=SIME?.EM  <03 

RHQ (5 , 1 3  -$CP  T  (PARA  fl  J  M  <?.  /pmo  (  Id  U  -(!«/££ £M(1))}  ) 

RHO  (6#I>  ®V1MC*A-CC*  A  !. •?-)*■ '}  t  -  (  -  A*CC*  ft  A  *  C  *  ° ) 

RH05?  ,D  *V!»  CC#AA*r.r»  ^*^>*Vt*  \  ->u*cc  *A*  c*  1 ) 

33  RHOU,  n  ®V1*CC’TT*  V2*c*30 
RETURN 
END 
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20 


SUBROUT  I NE  PCLY1.E  G  (II) 

COMMON/ OLK  A/RHCC  8,  60)  tELEM  (10)  ,  PARA(  10)  ,6(60) 

COMMON/ DLKB/F  (8,  5)  ,OP  (6,5)  ,DO.  (  6 ,5)  ,0  R  (8  ,  5 )  * RR  ( 6 ) 

A=RHO (4, II )/RhO( 1, xw 

6=A*A 

C=A*B 


Dr  A*C 
E=A*D 


*(i 

<1 

(1 

(i 

(i 


f~a*e 

AA=SORT { 1. -B) 

BB=1»-B 
CC=AA  *8B 
QD=BB*BB 
£E=8B*CC 
FF=BB*OD 

AAA=-RH0(4,II)*RH0  (2,11)  /RR (3) 

B8B=~  RHQ  (4,  II  )*RHO  (3,  II)  /RR (3) 

!CC=(RR(2)  -RFC  (4,  ID  *RHC(4  »II)  )  /RR(3) 

‘  ,i)  =  1.  5*E“Q.  5 
,2)=2,5*C-1.3*A 
,3)=0.125*  (35  .*0-3O.*B+3.) 
i 4)  =  G . 1 25  *  (63  .  *  E-7  U  .  *0  15,*  A  ) 
5)=0.0625*<231  .*F-3i5,  *D+i  0  5.*B-5.) 

1)  =  AA*3,*A 

2)  =AA*  (7 .5*1-1 *5) 

3  )  =  A  A  *  (17  .5*C-7.5*A> 

4) =AA*0,  1  25*1  (3  15,  *0-21  G.*-3+l5.  ) 

5)  =AA*G. 125*  (603.  *E-63  G.*C»i05. 
i>=88*3. 

2)  =  39*15.*A 

3)  =03*  (52*5*  3- 7*5) 

4) =BB*  (157.5*0-52, 5* A) 

5)  =  BB*0.  125*  (3  465,  *0-1  890.  *6*105.) 

2 )  =CC*15  * 

3) =CC*i:5.*A 

4) =CC* (472.5*3-52. 5) 

5)  00*  (  1732.  5*0-4  7  2.5*  A) 

3)  =00*  105  . 

4 )  =00*  94  5 • *A 

5 )  =00*  (5667,5*6-472.5) 

4>=EE*9A5. 

5  )  =EE*  10  395.  *  A 
5 )  =  FF*  10  295. 

I-  i»  7 


A) 


20 

20 

-1 


J=l#5 


■l 

J)=CCC*(  (P<L,  J)/AA> 
L,J)  =  OR  (I  ,J)*AA  A/CCC 
J)=OR  (I,  J)*83B/CCC 
'INUE 


•(  K*A*  P  (I,  J)  /BB)  ) 

Reproduced  from 
best  available  co 


copy. 
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SUBROUTINE  ?  CNAL  ( ID 

COHMON/BLKA/fiHO(3,  60)  ,EL£M  (10)  ,  PARA(  10)  ,E  (60) 
COMMON/ ELK  8/ P  (3, 5)  ,D» (6,5)  ,CQ( 8 ,5)  ,0  R(6  ,5)  ,RR(6> 
r'OMMQN/0LKC/BJ<5)  ,  C(  7, 5)  ,S  (7,5  )  ,FAC(  7, 5  ) 
COHHON/PLKC/FERT (7,3,60) 

PX=PARA(i)/RR(3) 

OOPX^Oo 

ODPY--0. 

DDPZ=  0* 

DO  20  1=1,5 
BB-BJ(I)/RR(  I  +  i) 

CC= (1  +  2) *P  (i  ,1) 

00=DP(1,I)  *RR(2) 

EE=OQ  (1, 1)  (2) 

FF=DR(i,I)  *RR(2) 

QOPX  =  DO*  (CC*  ':.HO(  2,  ii  ) -OO)  +  CDPX 
ODPY«DB*(CC*?hO(3,II)-EL)  +  DDPY 
DOPZsOB* (CC*RH0(4, II »-FF)+DDPZ 
20  CONTINUE 

PERTH  ,1,11)  =FX*ODPX 
PERT (1,2,11) =  PX*DDPY 
PERT (1,3,11) =  PX*  DO  PZ 
RETURN 
END 
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S UBROUTINE  Tr.SSER  <AL,II) 

C’OMMON/nlKA/RFO(  3,  60  )  ,ELEH  (10)  ,PARA(1Q)  ,E(6G) 
COHHON/BLKQ/F  (6,  5)  ,DP(8,5>  ,DQ(  8,5)  ,DR(8 ,5)  ,RR(  6) 
COHMOM/BLKC  / B J (5)  ,  C(  7, 5)  ,S  (7,5  )  ,FACt  7,5) 
COMMON/BLKC/FERT (7 ,3,60) 

PX=P AR A ( 1 )  /RR(3) 

RRR=RR  ti)*  (RHH2,  II)  *KH0(2 ,11)  *RH0(3 ,11)  *RHC(3,  ID) 
PXX=-PARA (i)  *RHQ  (3  ,11)  /RRR 
PYY=PARA(i)*KHO(2 , II ) /RRR 
ODPXl=  0 » 

OCPX2=0» 

DOPY  1=  0 • 

DOPY2=0. 

DDPZ=0. 

00  20  1=1,5 
K-I  +  l 

DO  20  J=2,K 
H=J- 1 

C01  =  C0o  (M*  AL  ) 

S Il-SI N( M*  AL ) 

CCSI1=C01*C<  J  ,  I)  +SIi*S  < J , I  ) 

C0SI2  =  C01  +  S(J,I)-SI1*C(J,I) 

AA~-<K  +  1)*P(J,I)/RR(  1  +  1) 

BO=DP(J,I)/RR(J.U) 

CC=OQ(J,I>/RR(I+l) 

OD*OR(J,I)/RR(Ifl) 

EE=M*F(J,I)/RR(I  +  1) 

ODPX 1=  ( A  A  *  RHO  ( 2,  II)+B3*RR<  2) )  *  COSI  If  OOP  X 1 
OOPY1=  (AA+fiHC  (3,  II  )  +  CC*RR(  2))*  COSI1+  DOPY1 
D0PZ= (AA+RHO  (4,11)  +00* RR  (2  )  )*C  OSIi+DOPZ 
Q0PX2  =  tE*C0SI2+QDPX2 
20  CONTINUE 

OOPXi  =  PX*OCPX'l 
OOPY2=PYY*OCPX2 
OOPX2=PXX*OOFX2 
DDPY1=PX*0QPY1 
PERT (2,1,11) =QCPX1+Q0PX2 
PERT  (2,2,11)  =C0PY1+00PY2 
PERT  (2,3,11)  =PX*  DQPI 
RETURN 
END 


ga/kc/73-? 


SUBROUTINE  SECTOR  (At,  II) 

COMMCN/BLK A/ FPO(  3,60)  »  t  L  EM  ( 1Q>  ,PARA(  10)  ,E<60) 
COMMON/BLKU/F  (8,  5)  , DP  (8,5)  ,  00. (  8  .  ^ )  , 0 R  (8  ,  5 )  ,  RR  (  6 ) 
COMMON/ OLKC/EJ (5) , C( 7, 5) ,S (7,5 ) ,F&C( 7,5) 
COMMON/PLKG/PERT (7,3,60) 

PX= PARA(l)  /RR  (3) 

RRR=RR  ( i  >  *  (RHO(2,II)  *RHC(2  ,11)  +RH0(3,II>  *RH0(3,II)  ) 
PXX--PARA (i) *  RPO (3,II)/RRR 
PY»=PARA(i)*RH0(2,  ID/RRR 
ODt  X1=0 ♦ 

ODPX2-0* 

COPY  1=0* 

OOPZ1=0. 

KKK=i 

00  20  1=1,5 
K-I  +  l 
J=I  +  2 

KKK=(2*K-i ) * KKK 
COl-COS ( K* AL ) 

SI1=SIN(K*AL) 

C0SI1  =  C01*C(  Jfl)  *SIi*  *,I) 

C0SI2  =  C0i*S  ( J,I>  -SI1*C  ,1  ) 

AA-“KKK^(K  +  l)'f'COSIi/RR(<> 

ODPXi=AA*RFC(X,II) +OCFX1 
DDPYi=AA*RHO(3,II) +DQPY1 
0DP71=AA*RP0  (4,11) +D0PZ1 
ODPX2-  KKK^K*  CCSI 2/  RR  ( i< )  +C0  PX2 
20  CONTINUE 

PERT  (3,1,11)." FX* CD ~n  I  +PXX* DOPX  2 

PERT (3,2,11 5  sFX*03PY1  +  PYY*00Pa 2 

PERT  (3,3,11)  sFX*'OOPZl 

RETURN 

END 


SURE  OUT  I NE  t  UJN  A \J  V 

C0MMCN/BL<A/«P0(d,6(n  » El  EM  (10)  ,PARA(  10)  ,E(6  0) 
COMMON/BLKC/PERT (7,3,60) 

COHHOn/BUKF/ CRS(  3 , 6G ) ,RRM( 3,6G 1 
PARL  =  0 «  Q  1 226  A'-l  A  (  1) 

A=RHOC2,I)-RRH(l»I) 

0-RHO ( 3 , 1 ) ®R  R  M2 , 1 ) 

CsRHMA,  I)  -RRM(J,I) 

0=SQRT<A*A+fc*P+0*C) 

0*0*0* 0 

00=S0RT(RRHU,I)  *RRM<  i»t>'  RRM(  2 ,1)  *RRM<  2,1) 

1  ♦PRM(3,1)*RPM(.?,I)) 

EE*DD*0D*DC 

PERT (4,1,1)  =  -PARL* (A/Q  +  RkH(1,I  ) /£E) 

P£RT(4,2, 1)  =  -PARL*  (P/D+RRH  (2,1  )  /EE) 

PERT  (4,3, 1)s-PARC*  (G/C+KRH  (3,1.  )/££) 

RETURN 

END 
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SUBROUTINE  SCLAR  (It 

COMMON/dLK  A/RHOC  3,  faC  )  ,  EL  EM  (13)  ,PARA(  10)  ,E(60> 
CCMMON/PLKC/FERT (7 ,3,60) 

COMMON/BLKF/RRS(3,  60  )  ,  RRM(  3 ,60  ) 

PARS=33295Z. *FA*A ( 1) 

A=RHO  (2  9 1 )  -RES  (1 , 1  ) 

B^RHO  (3,1)  -RRS  (2,1) 

C  =  RHO(A,I) -RRS<3,I  ) 

D=SQRT(A*A+8*SrC*C  ) 

D=D*D*0 

DD-SQRT  (RRS(  1,1)  *RPS  ( 1 , 1 )+ RRS  <  2  ,1 )  *RRS  ( 2 , 1 ) 

1  +RRS (3,1) *RRS  (3,1)) 

EE=0O*CD*OC 

PERT (5,1,1) s-PARS* (A/C  +  RRS  (1,1) /EE) 

PERT (5, 2, 1) a -PARS* (B  /  0  +  RRS  (2,1) /EE) 

PERT  (5,3,1  )  =  -RARS*  (C/O  +  RRS  (3,1)  /EE) 

RETURN 

END 


suoRoun e_ rac iat  t n 

COHHON/blKA/RhO(3,  6u  )  ,  E  L  EM  (10)  ,PARA(  10)  ,£(60) 
COMMON/ fllKO/FERT (7 ,3,6C 
COMMCN/BLKF/SRS(3,60) ,RKH( 3,60 ) 
AsRHO(2,I)-R«S(l,  I  ) 

B^RHO (3,1) -RRS(2,I  ) 

C^RHO (4,1) -RRS (3,1 ) 

RSOcSCkT (A^A  +  B^B  +  C^C) 

6AH*ASINU  ./RPC(1»  I)  ) 

£PS=ACOS  ( (  A*  RHO(  2  »  I)  ♦°*RHO  (3,1)  +C*RHO(4,I)  ) 

1  /<RSQ*RHG(1,I  )) ) 

AEPS=ABS(EFS) 

PAR*  (PAPA<  10)  *P ARA  (9)  )  /  (RSO*RSD*RSO) 

IF  (AEFS*G  E,  CAH)  GO  TO  50 
GO  TO  60 

50  PERT(6,  l,I)aPAR+A 
P£RT<!‘.2,I  )  =  FAR*B 
PERT (6|3,i*=FAfi*C 
RETURN 

60  PERT ( 6, 1 , 1 )  =  0  » 

PFRT (6,2,1 ) -0. 

PERT ( 6, 3, i ) -  0. 

RETURN 
END 
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SUBROUT INE  TAEL Ed  I ,  JJ,N,N0SC) 

COMMON/ BLKO/PERT (7 ,3. CO) 

COMMQN/BLKE/ [  IFF  £  i  0,  6  0  ) ,  TA  e  <4,  6  0)  ,  H<  3 , 5  )  ,TA3S  ( 2 , 7 , 3  ) 
M=N-i 

00  10  1=1,10 
DO  10  J= 1 , 5 
OIFF ( I ,  J )=  G  « 

K=N+i-J 

10  DIFF<I,K>=0. 

DO  20  1=1,5 
12=2*1 
11=12-1 
DO  18  J  = 1 , M 
IF  (I.GT.l)  GC  TO  17 

OIFF  ( 1 ,  J  )  =  PE  RT  (II,  JJ,  J  +  D-PERT  (II,  JJ,J) 

GO  TO  13 

17  OIFF  (Il,J)=QIFF(Il-l,J+i»-DIFF(Il-l,J) 

18  CONTINUE 

DO  19  K=2,N 

19  OIFF<I2,K>=OIFF(Il,K)-OIFF  (11,  K-l) 

•20  CONTINUE 

TA8  (4,0)  =TAOS  (2,  II,  JJ)  +PERT(II  ,  JJ,&)  *0,5 
TA8(3,6)=TA8S(1,  II  ,Jvi)  -PERT  (II  ,  JJ,6)  /12* 

DO  30  1=1,5 
12=2*1 
11=12-1 

TAB(3,6)=TAB  (3,6) -OIFF  (12,  6>*H  (2,1) 

30  TAS  (4,6)  =  TAR  (4,6)-  (OIFF  (II,  5)  +  OIFF  (I  1, 6)  )  •U.5  +  H(i,n 
DO  40  1=7,  N 

TAB (4 ,1}  =TAF  (4,1-15 +  PERT  (I  I,JJ  ,  I) 

40  TABU,  I!  =TAG  (3,1-1)  + TAB  (4,  I-i) 

00  50  1=1,5 
J=6-I 

TAG(4,J)  sTAE(m,J  +  1)-PERT(II,JJ  tvl+l) 

50  TAO<3,  J)=TAB  (3,J  +  t)-TAfi(4,  JJ 

A=DIFF  (10,7)  -CIFF<  jiQ  *  6> 
a=DIFF(10,N-55-OIFF;  OjN-6) 

DO  60  1=1,5 
J=6-I 

DIFF  ( 1  0,  J)  =OIFF(lO,J  +  i)-A 
J=N“5+I 

60  0  IFF ( 1 0 , J) =CIFF{  io ,  1)  +  8 

DO  70  I-i, 9 

jno-i 

L=J/2 

L L*N* (J~i)/2 
LLL=  ( J+ 1 )  / 2 
DO  71  U=Li,H 

r i  dtp f<  j,  t  1 ;  r  i"-;  j ,  ti - 1 )  +8T rr  ( '  *  1, 1 l/llu 

IF  (I,  to.  9)  GC  To  70 
00  72  11=1, L 
12=1  +  1-11 

7  2  DtFF  (J ,  I?)  =OIFF(J,  IS  +  i  )  -01  FF(  J  +  1, 1  2+  1-L/LLL) 

70  CONTINUE 

00  80  1  =  1, N 
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TflB(i,I)=TAP(3,n  *  PERT  (II,  JJ,I  }  /12. 

TABC2,  I>=TAP  (4,1) -PERT  III,  JJ,I)  *0.5 
00  80  J=i,5 

TAO(i,I)  =  TfiG  U, I)  +DIFF'2*J  ,I)*H  (2,  J) 

IF  { I  *  £  Q  •  1 )  GO  TO  60 

TAB(2fI)=TAE(?,I>  +  <91fF£2*J-l,  I )  +DIF  F  ( 2*  J-1  yi-i)> 
1  *0.5*H(lfJ) 

80  CONTINUE 

TABS<1,II,JJ>=TA=>{1,6) 

TA0S(?,II»JJ)=TA3{2,6) 

RETURN 

END 


SUBR PUTIN E  I  h PROVE  (  I»  N_i Hi  NN ,MO SC) 

COMMON/flLKA/ Shot  3,  GO)  ,’LLGH  ( I  J)  ,  PARA  £  iO)  ,E(60) 
CONHOH/BLKC/FERT (7 ,3 , GO) 

COMMCN/BUKE/CIrFUCtSO)  ,  TA  B  (4,  60)  ,H(3,5>  .TA9S(2»7,3> 
C OMMCN/BLKG/  '<EF(  3, SO)  ,  EREF  (60)  ,  ELREF  (10) 

00  50  XTER=1,NN 

00  20  K=1,N 

X~REF  £  2 ,  K) +  P~ R T £  I , i, K ) 

Y-REr £3,K) ♦P£RT( I , 2, <> 

Is RfF C4»K)*PEKT( !♦ 3, K) 

XX*PE?  T  ( 1 , 1 ,  K  J *  < RE  F  <  2 ,  N)  40 , 5*n  E  RT  ( I,  1 ,  K  )  ) 

VVepERT  (1 ,  2  ,  it )  *(  RE  F  £  3  ,  K)  .  5*°  c  F  T  £  I,  2 ,  *  )  ) 

22-PERT  (I,  3,i0*0£rm,<)  ♦3.5»PL‘PT(  It  3 ,  K  i  ) 

Q«  E  XX+  Y  Y  +  2  2)  7  W:F  (  U  <)  *REF  (i,<  )  > 

.7. 5*  3*  ->  -IB.  3  25  *0*0*  0+36.  G?5  *Q*Q*0*0 

1  -ifi7.68?S*C*C»0*C*C  *4C2 ,1375*0*0*  0*C*C 

2  *854,04844  *0 *0 *0 *0 *0 *Q *0 

PARsPARA  (1  ?7  £R£~  ( l  ,K)  *  R  £  F  (  1  ,K>  *  REF  (1  ,F>  ) 

&ORfti*PAR*<F  *C*'<"RfRT  <  I » lr  X) ) 

C0Rft2*PAS* (F*C*y-RFRT  (1,2, K) ) 

C08R3*PA .»*  (F*C*2*^tRT  (1,3,  *>> 

PERT  'X  *lfX)sP£Rt  (7 ,1  ,n)  ♦v.ORRl*  W*H 
PCRT a,  ?#<)*■  FERT  (7  j?  ,  k  )  ♦  CO RR2*  N*W 
PERT  U  ,  3  ,  K  )  a  PERT  (  7  , 3  ,  K  >  *C0  R f>3*  K  »  w 
REF ( 2 } K ) *  X 
REF  (3,K)aY 
REF  (i,,  K)  =2 
20  CCNTT SUE 

IF  (ITER.EO.NN)  GO  TO  50 
OO  40  J=l, 2 
T ABS ( 1 , 1 » J ) ~  C • 

TAU$(2,I,J)=0, 

GAEL  TABLE  (I,J,MfNOSC> 

00  30  n*i,n 

30  FCRT (I, J,K)=TAB( 1, K) 

40  CONTINUE 
50  CONTINUE 
RETURN 
EN'O 
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SUOR  OUT  I  HE  G  S_C  UL  AT  (NOSC  > 

COM H 0 N / i3 L <  A / •< P  0  (  -i ,  GO)  ,tLEM  (10)  ,PARA(10)  ,£(60) 
GOMMON/ULKC/FERT  £7  ,3,60) 

VX=RHO (6,NCSCi 
VY=RHO {  7 ,NOSC ) 

VZ~RHO(  B,  NOSC) 

V2=VX*  VX+VY*VY*7Z*  VZ 
X=RHO (2, NOSC ) 

Y=RHO (3, NOSC ) 

Z=RHO ( 4 , NO  SC ) 

R-SQRT (X  *X +Y  *Y+Z  *Z ) 

ELEM(i)  »1,  /  (2./R-V2/PARA  '1 3  ) 

HX=YJf-  VZ-  2*  V Y 
HYsZ*VX-X*  VZ 
HZ=X*  VY-Y*  VX 
H2-HX*  HX  +  H  Y*  FY*H  Z*  HZ 
£L£H ( 6) =H2/PARA< 1) 

£LEH(?)=30kT  <l.-cLEM<6)/£LEM(i)  ) 

A  =  V2"*PAR  A  (  i  )  /R 

B=X*VXfY*VY+7*VZ 

£X-  ( A  *  X-t?*  VX  )  /  PA  Q A ( 1 ) 

EY*<A*Y-6*VY)/PARA  (1) 

F7#(A»Z-PvVZ)/f,ASA  (1) 

AA^SOPT (H>> 

El.EM  (  3)  -  ACCS  (h?/  AA  ) 

ANsSQRT  (HX«H:<  +  HY*-(Y) 

£ttrt(  S>  -ACCS  (-HY/AU) 

IF  <HX.LT,  6*)  ELSH  (5)a2.*PARAC2)-C.L?-M(5) 

£LEM  {  0)  a  A  DCS  (  <*HY*£XfHX*£Y  ) /(A  N*El.LM  (?)  )  ) 

IF  (F2.U.0.)  EL£<(4)  =  2.*PARA<2)*ELEM<4) 

EL£H< A) a£LEM (4>*PACA (3) 

ELEM  <  5> -  EL  EH ( S ) *  PA  PA  < ? ) 

ElEM(3>sELEM(3^°APA(3» 

PARA<  4)=SQRT  ( f:  A-i  A(  j )  /  (  EL  EM  (li  *  CLEM  {  i  )  *ELGM  (1)  >  ) 
ECOSa  (EX*X  +  £Y»Y*-CZ*Z)  /  (FLc  m  (/.)  •  RJ 
IF  (AuS(£COS).l£.l.)  GO  TO  10 
PRINT  100,  LCCS 
100  FORMAT  t/,5X,Eli.6  ,/) 

IF  (ECUS  *GT«  0, )  ECOSsl* 

IF  .<ECOS.LT*  C* )  ECO^a-i. 

10  e (i) =acus<£ccs) 

IF  (O.LT.O.i  E  <i>  a?.  *PARft<  ?)-“  (1) 

£  (NOSC) -ACCS  { (a  EM  (?)  ♦  £CGS  )7(1  .  ♦LLEH  (2)  *€CC<) ) 

IF  (n.LT.O.)  t  (NOSC)  5  2.  *  F>  RA  (?  )  -£  (  NO  SC ) 

PARA  (  0 )  -  <1 ./  PARA  (  4 ) )  •  <E  ( NO  SO*  £  LEM  (2  )  *S  IN  (£  (HUSO  )  ) 

RETURN 

END 


90 


GA /mg/73-7 


S UBROU T I NE_ PE <? IG2E  (  N p FR _jU  ) 

COMMON/ OLKA/  RHO(  3 , 60  )  , El  EM  (10)  ,  PARA(  10)  ,G(60) 

U  =  0. 

A=RHO  ( 1  ,  NPER-2) 

B=RHO(i,NPER-i> 

C=RHO(i, NPER) 

0=RHO(i,NPER*l> 

£=RHO(ifKPER+2) 

A4=( A-4U*8*6. *C-4.  *D  +  E)/(2  4.) 

A3=(-A  +  2«*e-2,^rH-E)/{12.) 

A2=(-A+16.  *B-30.*C+16.*0-E  >/(24.) 

Ai  =  (A-8.*8  +  8.*C-E)  /  ( 1 2  • ) 

M=0 

tO  A4*ll*U*U  +  7*  *A  3*1.1*"  +  ?-  *A2*U+A1 

FF=12.*A4*U*U+o.  »A3*U  +  2.*A  2 
H  =  M+1 
U=U~F/FF 

IF  (M.GE.1Q)  GC  TO  20 
IF  (A8S(F)  .L2..033  00  0  1)  GO  TO  20 
GO  TO  10 
20  CONTINUE 

IF  (ARS(U)  .IT. .001  )  GO  TO  23 

IF  <  U  •  G  T  • ,  <3  9  9  •  AN  D «  U»  L  T  .  1,0  01)  GOTO  24 

if  <u.cr.-i.  cr.i.  v;r,  u.iTc- ,99?>  go  to  25 

AAA-UMU*U-i,)*{'J*U~H,)/24  . 

AA=AAA/ ( U  +  2 . ) 
er-*-<AAA*4,l  /  ( U+  1 «  ) 

CC= AA A*0 » /U 
00s*{AAA*4  .)/  (U-l.  > 

£E=AAA/ (U-2, > 

GO  TO  27 

23  DC=0, 

ec=o. 

CC=i, 

GO  TO  26 

24  CC=0. 

00=1. 

0B=  0  * 

GO  TO  26 

25 

cc=o. 

00=0. 

26  AA»0  ♦ 

EE=  0. 

27  CONTINUE 

00  30  j=i,e 

RHO(J,U  =  4A*«K><  J,  KPER-2M  FO*RHO(  J  , K  P£!>1 )  ♦ 
l  CC*PH0(J,NP£RH‘00*RHC(J,NPER*-  i)t£E*RHO  (J,NPER*2> 
30  CONTINUE 
RETURN 
END 
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SURRO UTTN F  r£ l_TAj/ 

COMMCN/3LKA/^hC(d,  £3  >  , EL  EM  (10)  ,  PARA(IO)  ,E(£U) 
COHMON/RLKG/«EF(  i,  00)  ,EREF  (60)  ,  ElREF  (10) 

RA=£REF  (19) 

VA=EF:EF(20  ) 

D3=£LEM(3)  -ELKEP'  (3  ) 

04=2«  *PARA  (2  )  ~EL  FM  (4) 

05-=  (ELREF  ( 9)  -EL- H  (  5)  )  *SIN{  EREF  (  33)  ) 

ONU=-ATAN?  (06  ,03  ) 

ENU=COS  (04  +  CM) 

VAZ=SQRT  ( (  F  AHA  (1 )  *  PA )  / (EPE  F (31 )  *RH 0(1,1)))* 

1  (l.fEFEF  (32)»ENU)/(1.  +  EREF  (32)) 

DVWsQ5*VAZ/SINONU> 

0 6=  EL  R t  F  ( 4  )  -  E  L. EM  (4  ) 

0  VN  s  0  6  v'  V  A  v‘  E  R  £  F  (2  1  '•  *F.Rf  F  {  22  ) 

DVN  =  DVN/(2  .  'EREF  (2  i)  *£RFF(  22)  -  RA) 

D7- ( 2«  *ELE  H(  «}-RHO(l  ,  1 ) f FR EF (4 1 ) ) /2, -EL E M ( 1 ) 

OVTi  =  PARA(l)  *07/  (  ?  «  *  E  R  E  F  (2  1>*FR£F<  21 )  *V  A) 

06=2.*  (EREF(  2  25-  1,  )*0\/Tl/V  A 
09-EL  REF (2) -EL Li (2  >-03 

RH  0(9,1)  =  SORT  (PARA  (1)  *  ( <  2.  /RHO  (  1, 1 ) )  -  (  1  .  /  (ELF.M  ( 1)  +0  7)  5 )  > 
0  y  T  2  -  R  H  0  ( 6  ,  1 )  ♦  0  -)  /  (  2 .  M  1 »  +  £  L  £  M  (  2  )  ♦  Q  3 )  ) 

01=  (2.  •  (FLER(  1)4-07)  *  (£L£h(  l)fD7  )*RHO  (5, 1  )*OVT2)  /PARA  (i) 

ElEM(iUtLEH<m01*0  7 

ECEM(2)“ELEM(2)*n?+P9 

ELEM (  3 )  »  1  Ilf. .w  (  ))  -02)  •  oiPAI  M 

ElEH(4)  =  (ELtM{4)  ♦06)  ’FAR A (  3 ) 

ElEM(5>s  (Ei£R(6)  OS)  *PARA{  3) 

6LEH(6)  a£LcM  (1)‘  (  1  .-EL  EH  (2  )  *£l  EM<2)) 

PAPA  (4)sS0RT  (PARA(  1)  /  (  £l£M  U)  *EL£«  <t )  *E  LEH  ( 1) )  ) 
C0N«637ai65. /360  0» 

£R£F (i)*OVW*CON 
EREF  <  3)  =  Ovn*CCm 
EPXF  (5)  =  OVU*Cr:N 
EREF ( ?) -DVT  Z *CCM 
EREF  (9)s0. 

£REF{2)aA0S(EREP(l))  +EPEFC  2) 

E PEF  (  4 )  ®  A  0  $  ( £  F  £F  (3  )>  ♦ES?£F(  4) 

EREF  ( C )  «  A”S  <  £F:iF  (5))  ttREF,  €> 

EREf  (  8 ) -b  A3S  t  £ 4 £F  C  7  )  )  «•  ERE  F(  8) 

00  10  1-1,4 

10  EREF  (  9)s  AfiS  (  EFEF  (2T-D»  *ER£f<  9) 

EREF  <  '  0)  «  EREF  (SI  ♦EREFflO) 

EREF  (11)  *-72.*VA*0VTl/S0RT  (PARA  (1)  *£RCF  (21) ) 

RHOU  ,i)5£l£l'(M  /(  l.*EU"«(2n 

RHO  (5,  1) -SCPT  (FA2A  (!)♦<(  2,  /RHO  ( 1,1  H  -U  ./ElEHil ) ) ) ) 
RETURN 
END 
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Glossary  of  jc  ;;  Variables ,  Tho  following  Glossary  is  a  list  of 
major  /arialues  used  in  the  computer  model  listed  in  the  previous 


section. 

Intermediate  vcariables  are  not  defined 

since  they  are 

functions 

of  the  major  variables. 

ALO 

Geographic  longitude  of  the  vehicle,  A 

C 

ALPHA 

Cumulative  right  ascension  of  the  vehicle 

AT 

Right  ascension  of  the  vc 

hide,  « 

BJ(1) 

Zonal  harmonic  coefficisr. 

ts,  Jn  (n  - 

1+1) 

c(.i.j) 

Earth  geopotential  harmon 

ic  coefficien. 

ts,  C 

'iik 

n  13  J+l 
k  -  1-1 

DATS 

Time  from  launch  (days) 

DXFP(I.J) 

Tabular  difforoi ices,  6P' 

n  «  1,  j  - 

i  1  —6 

I  .  11 
tJ'T 

for 

n-evei  ] 
r 

J  j 

. 

•  for  n-ol J  j 

»*(!..•) 

~x  F(V 0 

V  P(I,J) 

J 

) 

v.  P(l,j) 

1)U1 

Time  froa  perigee  to  apeyce  (hra) 

wz 

Tice  from  apogee:  to  perigee  (hrs) 

DVH 

Hartal  irsrulse, 

n 

DVT1 

Tangential  iapulce  at  aj<\ 

DVT?. 

ial  irpulso  at  perigee,  a v^s 

DVV 

Orthogonal  impulse ,  Av 

s(D 

Secsntric  anomaly,  E 

CbKK(X) 

Gesi-any^r  axis,  a  (>1) 

Argument  ci 

r  Jo 

rig so,  u 

Kceentr!  'Myt  0  (1*2) 

jRclir.it ion,  i  (1*3) 


Longitude  of  Ascending  , 

Kos*,u  U-5) 

Gossi-laUa.-  rectus,  p  (1«6) 


KBrK(I) 


Utility  array  ( eccentric  anomaly  of  roforer.oe  orbit,  orbital 
elennntu,  ole) 
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ELR£F(J>  Nominal  orbital  parar.-terr,  (initial  values  of  L*LS«(l)) 
*  t>eries  expansion  of  parameter  q  in  cJnckos'  Method,  f 


x> 

ISECT 


Joofficicncs  of  central  difference  into' rati  on  formulas 
Given  by  i'qc  (i'G )  oiid  ( jM j 

M  !  2J 

H(1.J)  *  (*,  ).  H(24J)  *  6. 

•*»  JL 


1SECT  Integer  re  rm  t outing  the  section  of  the  orbit 

I.TC7C  Integer  representing  the  cay  of  orbital  life  from  launch 

*  Numter  of  data  points  considered  in  each  integration  section 

Integer  rerresotu  in.g  data  point  nearest  to  the  apogee  pcin. 
NDAY  Total  number  of  days  cf  orbit  life  to  bo  rnalysed 

iu<  number  oi  iterations  in  backs  inurovuisent 

NOSC  Integer  represent irr;  pcint  of  orbit  osculation 

Ni'fil  Integer  representing  data  point  nearest  to  U*  psrigee  point 

r(I,J)  it  go  nine  pclynonis.lft,  F'  f  *  f  ^ 

K  *■  J  -  • 

I'AHA(l)  ( 2**1  ^  •’  ufh  gravitational  laranalor ,  ;»;« 

{1*2}  constant,  rUo“35G97'‘>) 

03)  ArUl#*  tic  consult ,  x$?.WrjV5 lyl 'deg, rad ) 

ONO  (N0n/)5 

(1*3)  Bdfferw.v-  :^.ar  distance  O.ft??>06?  sin) 

Co*  "toiit  \2y'‘y***Vu7  K  »**u, ) 

(  Ti-7).  ,-,  r  c 

U**1)  fi-a  at  ^cuiati-tn  joint,  frv.r  the  pun-ions  j*rl*or«  pel 
U*p)  Arra~M>  bass  ratio  of  vehicle,  A/r.  ^  r, 

( 1  **  ■*  in  v  -  os  ens  nro  e  t cr  \  y  *  f  XX p  j  o  rv'/nr  ye  *  null  A 

C  *  *’ 


*  ii't'  /  ■?  »  »  v  ’<  *  .  v  . 

'**  *  ^  V  «'  /  ♦’  £  *  *• 

v-Vdt  V-*  acceierri 

AV.-W  i.£  Cr  A*'.- 

(1*1} 

/■  C  .  !■'_'■  J  ^  V'  x  pA 

U~2; 

Vos--  -lal 

(J-l)  * 

0-3  5 

GeciC'-'ial 

02)  y 

(3-h 

lunar 

l>3)  -• 

}  *  -*;  7 

f.  0.5  ,w 

1  i  *  a  - 

(j>l  -  f 

(1-7) 

only 

F«trar.eter  q  in  ~r  -ho  -ethed 
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BEF(I,J) 

IiH(K) 

RHO(I.J) 


RM(I,J) 

RRM(l,j) 

BRS(I,J) 

RS(l,j) 

S(I,J) 

T 

TAB(J,I) 


TABS(l,I,j) 
TABS ( 2 , 1 ,  J ) 
THETAO 
TT 

u 


USAVE 


W 


Utility  array  for  temporary  storage  or  data  from  RH0(l,j)  and 
for-  storage  of  nominal  orbit  position/velocity  components 

Integral  powers  of  radial  distance,  rn  (n=K) 

Reference  orbit  position  and  velocity  components 
(1=2)  Px  (1=6)  px  (1-1)  P 

(1=3)  Pv  (1=7)  Pv  (1=5)  P 

(l~4)  (1=8)  Pz  J  **  orbital  point 

Daily  0  values  of  lunar  ephemeris,  (Jc'day) 

(lni)  Right  ascension,  a 
.  (I- •■2)  Declination,  6ro 
(I----3)  Semi-diamoter ,  s 

Position  components  of  the  Moon  for  the  orbital  point  J 

(■  )  \  (1-2)  4  (1-3)  V 

Position  components  of  the  Sun  for  the  orbital  point  J 

(Itul)  *s  (1=2)  ys  (3>3)  as 

Daily  0  values  of  solar  ephemeris,  (d-day) 

(I“l)  x6  (1-2)  ys  (1-3) 

Earth  gcopotential  harmonic  coefficients,  S  ^ 

Time  (hr) 

Tabular  difforencos  and  sums  in  integration  table 


(J-l) 

6x^  (i  « 

1-6) 

C\i 

w 

fl 

*“3 

(i  “  1-6) 

(J“2) 

w6x^  (i  « 

1-6) 

(j»4)  z\ 

(i=I-“) 

n  “  J+l 
k  -  1-1 


Initial  values  for  6 r  (J  «*■  1,2,3  for  x,y,z~coinponente)" 

Initial  values  for  6ro  (l  “  1  as  in  PEKT(I,J,K)) 

Greenwich  sidereal  tine  at  vehicle  launch,  0 

’  go 

Time  (days) 

Total  time  spent  in  orbit  or  normalized  timo,  T 

Total  time  spent  in  orbit  prior  to  orbit  being  analysed,  (brs) 

Timo  interval  between  data  points  (hre) 
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